3
$\begingroup$

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:

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ Over a finite domain {x, a, b} or over all real numbers x? $\endgroup$ – Michael E2 Jan 21 '15 at 21:43
  • 2
    $\begingroup$ 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? $\endgroup$ – Michael E2 Jan 21 '15 at 21:45
  • $\begingroup$ 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 $\endgroup$ – Remi.b Jan 21 '15 at 21:50
3
$\begingroup$

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.

enter image description here

$\endgroup$
4
$\begingroup$

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

From Why `FindMaximum` doesn't work in my example:

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}}
*)

Mathematica graphics

$\endgroup$
3
$\begingroup$

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}]

enter image description here

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

$\endgroup$
3
$\begingroup$

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]]

enter image description here

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

enter image description here

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.

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]}]]

enter image description here

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]

enter image description here

... 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. *)
NumberLinePlot[fd2, {x,-3,3}, ImageSize->500, ImagePadding->20,
     PlotStyle->Directive[{Red,Thickness[.01]}]]

enter image description here

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]

enter image description here

$\endgroup$
2
$\begingroup$

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]]

enter image description here

$\endgroup$
0
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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