How to represent in a 3D graph the set of solutions to this system of inequalities?

Question

I am studying this piecewise function:

Piecewise[{{f1[x, h], x < a}, {f2[x, a, h], a <= x < b},     
{f3[x, a, b, h], b <= x < c}, {f4[x, a, b, c, h],  x >= c}}]

where

f1[x_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 cdf1[x])/(h + (1 - h) 2 cdf[x])
f2[x_,a_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (1 - cdf1[x] + cdf1[a]))/
(h + (1 - h) 2 (1 - cdf[x] + cdf[a]))
f3[x_,a_,b_,h_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (cdf1[x] - cdf1[b] + cdf1[a]))/
(h + (1 - h) 2 (cdf[x] - cdf[b] + cdf[a]))
f4[x_,a_,b_,c_] := pdf1[x]/pdf[x] (h + (1 - h) 2 (1 - cdf1[x] + cdf1[c] - cdf1[b] + 
cdf1[a]))/(h + (1 - h) 2 (1 - cdf[x] + cdf[c] - cdf[b] + cdf[a]))

where, in turn

pdf1[x_] := PDF[NormalDistribution[1, 1], x];
pdf[x_] := PDF[NormalDistribution[0, 1], x];
cdf1[x_] := CDF[NormalDistribution[1, 1], x];
cdf[x_] := CDF[NormalDistribution[0, 1], x];

I want the first and the third piece of the piecewise function to be smaller than a certain value $l$, say $l=0.48$; and the second and fourth piece of the function to be greater than the same $l$. The following code represents the intervals of $a$, $b$ and $c$ such that the four inequalities are satisfied:

{amin = NArgMin[{f2[a], f2[a] > l}, a], amax = NArgMax[{f1[a], f1[a] < l}, a], 
a = RandomReal[{amin, amax}]}
{bmin = a, bmax = NArgMax[{f3[b], f3[b] < l}, b], b = RandomReal[{bmin, bmax}]}
{cmin = b, cmax = NArgMax[{f3[c], f3[c] < l}, c], c = RandomReal[{cmin, cmax}]}

Note however that, for example, you have to give a random value to $a$ since the lower bound of the interval of the $b$'s depends on $a$.

I would like to plot in a 3D Graph all the possible combinations $(a,b,c)$ such that the system of inequalities is satisfied. Moreover, I would like to implement the Manipulate command so that I can study how the set of all possible combinations $(a,b,c)$ changes as $h$ varies between 0 and 1. If a 3D graph was not feasible, actually any graphical representation of the set of solutions of $(a,b,c)$ would be fine anyway. Thank you if you can help.

Answer 1

The requested plot can be obtained as follows. First, compute the limits on a, using the expressions in the question.

{amin = NArgMin[{f2[a], f2[a] > l}, a], amax = NArgMax[{f1[a], f1[a] < l}, a]}
(* {-0.233969, 0.544866} *)

Next, compute the corresponding limits in b.

btab = Interpolation@Table[{a, NArgMax[{f3[b], f3[b] < l}, b]}, {a, amin, amax, .02}];

and finally in c.

ctab = Interpolation[Flatten[Table[{{a, b}, Quiet@NArgMax[{f3[c], f3[c] < l, b < c},
    c]}, {a, amin, amax, .02}, {b, a, btab[a], .02}], 1], InterpolationOrder -> 1];

Since a < b < c, the plot is given by

Plot3D[{b, ctab[a, b]}, {a, amin, amax}, {b, a, btab[a]}, AxesLabel -> {a, b, c}, 
    LabelStyle -> {Black, Bold, Medium}, BoxRatios -> Automatic, ImageSize -> Large, 
    Mesh -> None, PlotStyle -> Opacity[.5]]

where the allowed values of c lie between the upper and lower curves (both of which have the same projection on the {a, b} plane). The limits on a, b, and c were obtained here.

Addendum

As noted by kgir in a comment above, RegionPlot3D also can be used once the logical structure and interpolation functions have been determined. RegionPlot3D is much slower here than Plot3D but produces a figure that some may prefer. The choice is a matter of taste.

Quiet@RegionPlot3D[amin < a < amax && a < b < btab[a] && b < c < ctab[a, b], 
    {a, amin, amax}, {b, amin, btab[amin]}, {c, amin, ctab[amin, btab[amin]]}, 
    AxesLabel -> {a, b, c}, LabelStyle -> {Black, Bold, Medium}, 
    BoxRatios -> Automatic, ImageSize -> Large, PlotPoints -> 500, Mesh -> None]

Second Addendum

As noted by the OP in a comment below, the computation of ctab given here often fails with the error

Interpolation::femimq The element mesh has insufficient quality of 0.`.

as h is increased. Basically, the triangular grid of Interpolation becomes tangled near the edge at which the two curves in the first plot come together. A simply way of eliminating the error is to replace 0.02 by 0.06 in the computation of ctab. However, doing so decreases the accuracy of the plot significantly for 1/3 < h < 1. A much more accurate but slower computation involves replacing the three lines of code immediately preceding the first plot by

btab = Interpolation@Table[{a, NArgMax[{f3[b], f3[b] < l}, b]}, {a, amin, amax}, 
    (amax - amin)/20}];
fc[a0_, b0_] := Quiet@NArgMax[{f3[c], f3[c] < l, b < c} /. {a -> a0, b -> b0}, c]
Plot3D[{b, fc[a, b]}, {a, amin, amax}, {b, a, btab[a]}, AxesLabel -> {a, b, c}, 
    LabelStyle -> {Black, Bold, Medium}, BoxRatios -> Automatic, ImageSize -> Large, 
    Mesh -> None, PlotStyle -> Opacity[.5], PlotPoints -> 10, MaxRecursion -> 0]

The resulting plot for h = 1/2 is