What is the range of values $x$ for which $f(x)$ is higher than $k$ over a given domain?

What is the easiest way to answer the following question in Mathematica:

Given a function $f(x)=y$, what is the range of values $x$ for which $y$ is higher than some number $k$ over the domain of $x$ $[a;b]$?

For example consider the function:

What is the range of values of $x$, for which $f(x)$ is greater than $k=4.8$ over the domain [-2;2]?

As I am aiming to apply this method to a complicated function, I would like numerical approximations rather than an analytical solution.

• Over a finite domain {x, a, b} or over all real numbers x? – Michael E2 Jan 21 '15 at 21:43
• Note Reduce[-x^3 + 3 x^2 - 2 x + 5 > 48/10, x] works on the sample function...but I'm not sure what kind of difficult functions you have in mind. Would numerical approximations be satisfactory, or do you desire symbolic solutions? – Michael E2 Jan 21 '15 at 21:45
• Oh good points. Well, I will apply this method to a finite domain (of real numbers only). So let's say over the domain [-2;2]. And I would need numerical approximations. I edited my post. For info, my function of interest looks more or less like that Log[(E^(4 n qhat s) (1 - qhat)^(-1 + 4 n \[Mu]) qhat^(-1 + 4 n \[Nu]))/ NIntegrate[ E^(4 n qhat s) (1 - qhat)^(-1 + 4 n \[Mu]) qhat^(-1 + 4 n \[Nu]), {qhat, 1/(4 n + 1), 1 - 1/(4 n + 1)}, MaxRecursion -> 12]]. Thanks – Remi.b Jan 21 '15 at 21:50

Let

f[x_] := Sin[3 x^2] + x;


and criterion k = 1, all in the range [-2,2].

One can use Solve[], NSolve[] or FindRoot[] to get crossing points, but it helps to search for solutions using a number of starting points.

myCrossings = Sort@DeleteDuplicates[(x /. FindRoot[f[x] == 1, {{x, Range[-2, 2, .4]}}])]


Then, include the lower and upper limit of your x range:

PrependTo[myCrossings, -2];
AppendTo[myCrossings, 2]


Then, depending upon whether the value of f[x] at the lower limit of the range is greater or less than the criterion (here, k = 1), choose pairs of crossing points for the corresponding ranges of x:

If[f[-2] > 1, Partition[myCrossings, 2, 2], Partition[Drop[myCrossings, 1], 2, 2]]


(* {{0.443523, 1.02785}, {1.39912, 1.86883}, {1.94362, 2}} *)

If you want a more easily readable output:

#[[1]] < x < #[[2]] & /@ If[f[-2] > 1, Partition[myCrossings, 2, 2], Partition[Drop[myCrossings, 1], 2, 2]]


(* {0.443523 < x < 1.02785, 1.39912 < x < 1.86883, 1.94362 < x < 2} *)

It is true, in the extremely unlikely case that you have a local maximum or local minimum of f[x] at a root--and thus f[x] remains below the criterion or remains above the criterion on both sides of the root--you might have to check for solution ranges.

The easiest way most likely depends on the nature of the function. The function referred to in a comment by the OP seems to be like one in another post that it is concave down. The solution in such a case is quite simple and would consist at most of one interval. However, if the interval is to be $[-2,2]$, then it is not exactly that function. An oscillatory function is quite different, and the easy solution to a concave/convex function would probably not solve the oscillatory case.

Dealing with the more general case, then, I can think of two "easiest" solutions. One involving NDSolve is likely to be quite robust. The other uses Plot, may be less robust, but is likely to be faster. For highly oscillatory functions, some tweaking of the number of PlotPoints might be required. (An extremely oscillatory function is likely to require some manual adjustment in all cases, anyway.) Both methods assume that the function is continuous over the interval. A discontinuous function may be handled by applying the methods here to the intervals over which the function is continuous.

NDSolve method

higher[f_, k_, {a_, b_}] :=
Module[{ifun, in, s, t, x1},
Reap[
NDSolve[{ifun[t] == f[t],
in[a] == Boole[f[a] > k], x1[a] == a,
s'[t] == 1, s[a] == a,
WhenEvent[f[t] > k, {in[t] -> 1, x1[t] -> t}],
WhenEvent[f[t] < k, If[in[t] == 1, Sow[{x1[t], t}, Interval]]; {in[t] -> 0}]},
ifun, {t, a, b}, DiscreteVariables -> {in, x1}],
Interval][[2, 1]]
]


[Note: I may have found a bug using Reap with NDSolve. I had to use a tag with Sow/Reap here, because extra data was sown. I chose Interval for the tag; anything would have done.]

Plot/zero-crossings method

I used Mr.Wizard's adaptation of David Carraher's zero-crossings function, which can be found here: Find zero crossing in a list.

davidZC2[l_] := SparseArray[#]["AdjacencyLists"] & /.
SApos_ :> With[{c = SApos[l]}, {c[[#]], c[[# + 1]]}\[Transpose] &@
SApos@Differences@Sign@l[[c]]]

higher2[f_, k_, {a_, b_}] :=
(plot = Plot[f[s], {s, a, b}, PlotRange -> All];
With[{xy = First@Cases[plot, Line[p_] :> p, Infinity]},
Partition[
Join[
If[f[a] > k, {a}, {}],
Block[{s},
s /. FindRoot[
f[s] == k,
{s, Mean[xy[[#, 1]]], Sequence @@ xy[[#, 1]]}
] & /@ davidZC2[xy[[All, 2]] - k]
],
If[f[b] > k, {b}, {}]
],
2]])


(The plot is stored in a global variable for inspection in case the result seems questionable.)

Examples

n = 200000;
k = 2 n - 50;
μ = 10^-5;
ν = 10^-5;
qhat = k/(2 n);

Clear[fun];
fun[s_?NumericQ] := Log[(E^(4 n qhat s) (1 - qhat)^(-1 + 4 n μ) qhat^(-1 + 4 n ν))/
NIntegrate[E^(4 n qhat s) (1 - qhat)^(-1 + 4 n μ) qhat^(-1 + 4 n ν),
{qhat, 1/(4 n + 1), 1 - 1/(4 n + 1)},
MaxRecursion -> 12]]

higher[fun, 0.1, {0, 2}]
(*  {{0.0109522, 0.266125}}  *)

higher2[fun, 0.1, {0, 2}]
(*  {{0.0109522, 0.266125}}  *)


Random example:

g[x_?NumericQ] := NIntegrate[x Sin[20 t] Exp[-t^3], {t, 0, x}];
higher2[g, 0.05, {-2, 2}]
Show[plot, PlotRange -> {{-2, 2}, {0.05, .1}}, Frame -> True]
(*
{{-1.93866, -1.78522}, {-1.6299, -1.48044}, {-1.31331, -1.18644},
{-0.963237, -0.923499}, {0.738644, 0.853318}, {1.02542, 2}}
*)


One approach is to build a Piecewise function. For example, the piecewise function g[x] returns 1 whenever f[x] is in the region larger than c and 0 otherwise. Placed in a Manipulate, it shows the region over which f[x] is larger than the value of c.

f[x_] := Sin[3 x^2] + x;
g[x_, c_] := Piecewise[{{1, c < f[x]}, {0, True}}]
Manipulate[Plot[g[x, c], {x, -2, 2}], {c, -2, 2}]


The above example borrows the function f[x] from David Stork's answer.

Update: You can also use the options Mesh and MeshShading using the information from FunctionDomain[f]:

mesh = {Cases[fd, _?NumericQ, {0, Infinity}]};

Plot[f[x], {x, -3, 3}, PlotStyle -> Thick,
Mesh -> mesh, MeshShading -> {Blue, Red}, MeshStyle ->Opacity[0]]


For the second example below, using f2 and mesh2 = {Cases[fd2, _?NumericQ, {0, Infinity}]}; gives

Update 2: You can more conveniently use the result of FunctionDomain in setting the color function, e.g., ColorFunction -> Function[{x, y}, If[fd, Red, Blue]], in both examples in the original answer.

FunctionDomain (in versions 10+)

There are two new Version 10 functions, FunctionDomain and NumberLinePlot, that you might find useful.

Using the example in David's answer:

f[x_] := Sin[3 x^2] + x

fd = FunctionDomain[{f[x],-2<x<2 && f[x]>1.}, x,Reals]
(*  0.443523 < x < 1.02785 || 1.39912 < x < 1.86883 || 1.94362 < x < 2. *)

NumberLinePlot[fd, {x, 0, 3}, ImageSize->500, ImagePadding->20,
PlotStyle->Directive[{Red,Thickness[.01]}]]


Plot[f[x], {x, -3, 3},  Background->White,
ColorFunction -> Function[{x, y}, If[-2 < x < 2 && y > 1.,Red, Blue]],
ColorFunctionScaling -> False, PlotStyle -> Thick,ImageSize->400]


... and using OP's example:

f2[x_]:= -x^3 +3 x^2 -2 x + 5
fd2=FunctionDomain[{f2[x],-2<x<2 && f2[x]>4.8}, x,Reals]
(* -2. < x < 0.121115 || 0.790851 < x < 2. *)
PlotStyle->Directive[{Red,Thickness[.01]}]]


Plot[f2[x], {x, -3, 3},  Background->White, PlotRange->{0,40},
ColorFunction -> Function[{x, y}, If[-2 < x < 2 && y > 4.8,Red, Blue]],
ColorFunctionScaling -> False, PlotStyle -> Thick,ImageSize->400]


RegionFunction can also be used. With f[x] defined as in the answer by David Stork,

f[x_] := Sin[3 x^2] + x;
Plot[f[x], {x, -2, 2}, RegionFunction -> Function[{x, y}, 1 < f[x] < Infinity]]


You solve f(x)=k for x to get the borders of domains and then take a testpoint from each domain, eg. evaluate $f( (x_n+x_{n+1})/2)$ in order to see the type of domain. For solving have look at NSolve.

• I don't think NSolve will work on a function defined in terms of NIntegrate. – Michael E2 Jan 21 '15 at 22:36