# Why does the following Solve take so long?

Why does the following take Mathematica so long to reply to?

Solve[BetaRegularized[0.5, k, n - k - 1] == 0.035, k]


where $n = 1000$

By the way, In case there is no exact solution, how can I request Mathematica to find the nearest solution?

• I doubt you get an analytic solution. You can try with using NSolve. – José Antonio Díaz Navas Nov 25 '17 at 15:13
• For things like this, I think Solve should issue a warning something like: "Solve is best at finding symbolic solutions for exact polynomial equations. For this problem, FindRoot may be more appropriate if Solve cannot find a solution." – John Doty Nov 25 '17 at 15:32
• @JohnDoty Yeah, that would get us rid of many questions posted here ;o) – Henrik Schumacher Nov 25 '17 at 15:35
• @JohnDoty Note that Solve can handle several kinds of transcendental equations these days. Look at the docs, under "Scope", both "Complex equations in one variable" and "Real equations in one variable", for examples. – Michael E2 Nov 25 '17 at 19:36
• @Michael E2 sure, Solve has heuristics that can attempt to handle a variety of things, but its reliable core is exact symbolic algebra. – John Doty Nov 25 '17 at 23:20

Have a look for FindRoot. Per default, it tries to use Newton's method. However, the derivatives behave rather badly, so a secant method might be a more robust choice. Brent's method (mixture of bisection and secant method) works very well, here:

n = 1000;
b = 0.035;
sol = FindRoot[
BetaRegularized[0.5, x, n - x - 1] == b, {x, 0.4 n, 0.6 n},
Method -> "Brent"]
F[x] - b /. sol


{x -> 528.115}

-3.84415*10^-15

Note that 0.4 n and 0.6 n are supposed to bracket the root of the equation.

Solve originally was principally a solver of polynomial equations $p(x)=0$ and some equations like $p(f(x))=0$ for suitable transcendental functions $f$. But beginning around V7, I think, efforts have been made to extend the reach of Solve so that it might live up to its name. Many transcendental equations can be solved if the roots can be isolated, and if in no other way, their roots can be represented by exact, symbolic Root objects.

It often helps to have exact coefficients.

Here Solve works, if the root is isolated and a sufficiently high, but finite, working precision is used:

Block[{n = 1000},
Solve[BetaRegularized[1/2, k, 1000 - k - 1] == 7/200 &&
2 n/5 <= k <= 3 n/5, k, Reals,
WorkingPrecision -> 400]
]
(*  {{k -> 528.1147041429262349889...}}  *)


Alternatively, one can construct the corresponding Root object:

(sol = Block[{n = 1000},
With[{f = BetaRegularized[1/2, k, n - k - 1] - 7/200 & /. k -> #},
Root[{f,
k /. FindRoot[f[k], {k, n/2, 2 n/5, 3 n/5},
WorkingPrecision -> 400]}]
]]) // AbsoluteTiming
(*
{20.7221,
Root[{-(7/200) + BetaRegularized[1/2, #1, 999 - #1] &,
528.11470414292623498897865816253955484889... + 0.*10^-403 I}]}
*)


Remarks on the numerics of the equation

We can see that at machine precision, the function is constant for n = 1000 until the middle of its domain 0 < k < 999:

Block[{n = 1000},
With[{f =
BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # //
Rationalize},
InputForm /@ Table[f[N@k], {k, 0, 1000, 100}]
]]
(*
{0.965, 0.965, 0.965, 0.965,
0.9649999998691152, 0.4523748657765868, -0.03499999991364527,
-0.035, -0.035, -0.035, -0.035}
*)


A locally constant function generally stumps solvers that use local methods. One should consider RootSearch or findAllRoots; or perhaps SystemTRootsDumpGuessRoots (or SystemTRootsDumpGuessRealRoots) from this related Q&A: Faster variant of Reduce for finding zeros of holomorphic functions in a region

Perhaps it is not surprising that the derivative is difficult to evaluate at machine precision:

Block[{n = 1000},
With[{f = BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # // Rationalize},
Table[f'[N@k], {k, Join[{1}, 100 Range[9], {998}]}]
]]
(*
{-4.46485*10^155, -7.43223*10^246, -8.32217*10^275, 1.77529*10^275,
4.69065*10^254, -2.14015*10^216, 3.1861*10^253, -1.13516*10^276,
2.6163*10^275, 2.96348*10^246, -4.46485*10^155}
*)


You need ~500 digits to get accurate values throughout the whole domain; only ~400 are needed for k between 400 and 600, which perhaps explains why WorkingPrecision -> 400 works above.

Block[{n = 1000},
With[{f = BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # // Rationalize},
N@Table[
f'[SetPrecision[k, 500]], {k, Join[{1}, 100 Range[9], {998}]}]
]]
(*
{-2.79399*10^-300, -5.31457*10^-162, -7.31643*10^-86,
-6.31716*10^-38, -5.41614*10^-11, -0.0252418, -3.60875*10^-11,
-2.70387*10^-38, -1.82439*10^-86, -5.86766*10^-163, -2.79399*10^-300}
*)


I think there is something strange about the derivative of BetaRegularized.

Show[
Plot[BetaRegularized[0.5, k, 1000 - k - 1], {k, 450, 550},
PlotStyle -> Black],
ListLinePlot[{{450, 0.035}, {550, 0.035}}, PlotStyle -> Red]
]


I defined a function

f[k_] := BetaRegularized[0.5, k, 1000 - k - 1] - 0.035


and attempted to locate k where f[k] was zero.

I tried some numerical methods (e.g, FindMinimum, FindRoot) and had no luck.

I computed the derivative symbolically and numerically

D[f[k], k] /. k -> 520
(* -2.34977*10^206 *)


the result is inconsistent with the plot.

I defined an approximate numerical derivative

df[k_] := (f[k + 10^-6] - f[k])/10^-6

df[520]
(* -0.0108823 *)


I then used FixedPoint and a simple Newton's method to get an approximation of k.

FixedPoint[# - f[#]/df[#] &, 530, 20]
(* 528.115 *)


Here is a plot on a magnified scale near the solution

Show[
Plot[BetaRegularized[0.5, k, 1000 - k - 1], {k, 527, 529},
PlotStyle -> Black],
ListLinePlot[{{527, 0.035}, {529, 0.035}}, PlotStyle -> Red],
Graphics[{
Blue,
Opacity[0.6],
PointSize -> 0.025,
Point[{{528.115, 0.035}}]
}]
]


So my guess is that there is something strange going on with the derivative.

Note that I am completely unfamiliar with BetaRegularized so you may very well get a better answer.

• Very good point. When deriving BetaRegularized[1/2, k, 1000 - k - 1] with respect to k`, we see that fractions of pretty large numbers have to be juggled with, so that there is few hope to compute that accurately with machine numbers. – Henrik Schumacher Nov 25 '17 at 17:24