# Problem with NSolve

I would like to plot the solutions of an equation, for different values of a parameter. This is my code

cdf[x_] := CDF[NormalDistribution[0, 1], x];
pdf[x_] := cdf'[x];
cdf2[x_] = cdf[x]^20;
pdf2[x_] = cdf2'[x];

mix[x_, h_, p_] = h pdf[x - p] + (1 - h) pdf2[x - p];
ratio[x_, h_] = mix[x, h, 1]/mix[x, h, 0];
sol = NSolve[ratio[x, h] == 0.5 && (-5 < x < 5), x]
Plot[x /. sol, {h, 0, 1}]


For $$h$$ between 0 and 0.6 there should be three solutions, for $$h$$ above 0.6 just one. The plot, which runs in 5 minutes, returns only the solution for $$h$$ bigger than 0.6, while I believe I should get 3 different lines.

Can you help me fix this?

Thank you

cdf[x_] := CDF[NormalDistribution[], x];
pdf[x_] := cdf'[x];
cdf2[x_] = cdf[x]^20;
pdf2[x_] = cdf2'[x];

mix[x_, h_, p_] = h pdf[x - p] + (1 - h) pdf2[x - p] // Simplify;

ratio[x_, h_] = mix[x, h, 1]/mix[x, h, 0] // Simplify;

Clear[sol]

sol[h_?NumericQ] := NSolve[ratio[x, h] == 1/2 && (-5 < x < 5), x]


Generate data for a ListLinePlot. This is slow

hValues =
Join[Range[0, 0.1, 0.01], Range[0.15, 0.55, 0.05],
Range[0.551, 0.6, 0.001], Range[0.6, 1, 0.1]];

data = Table[{h, #} & /@ (x /. sol[h]), {h, hValues}];

data2 = GatherBy[data, Length];

ListLinePlot[{
Rest@Flatten[data2[[1]], 1],
Sequence @@ Transpose[data2[[2]]]},
PlotLegends -> Placed[Automatic, {0.75, 0.45}]]


EDIT: If you want to plot for various powers, parametrize the functions with the power. Further, as stated in the comments below, "the better solution is to use ContourPlot as suggested by @UlrichNeumann."

Clear["Global*"]

cdf[x_, pwr : _Integer?Positive : 1] := CDF[NormalDistribution[], x]^pwr;

pdf[x_, pwr : _Integer?Positive : 1] := (D[cdf[z, pwr], z] /. z -> x);

mix[x_, h_, p_, pwr : _Integer?Positive : 1] :=
h pdf[x - p] + (1 - h) pdf[x - p, pwr];

ratio[x_, h_, pwr : _Integer?Positive : 1] :=
mix[x, h, 1, pwr]/mix[x, h, 0, pwr];

ContourPlot[Evaluate[Simplify[ratio[x, h, #] == 1/2 & /@ {20, 10, 5}]],
{h, 0, 1}, {x, -1/4, 9/4},
FrameLabel -> (Style[#, 12, Bold] & /@ {"h", "x"}),
ImageSize -> 324,
MaxRecursion -> 5,
PlotLegends ->
Placed[StringForm["Power = ", #] & /@ {20, 10, 5}, {0.75, 0.5}]]


• Thank you for your answer, it's just great and very helpful. I've got just one question: why is there a little 'hole' between line 3 and 4? I would like to fix this if it were possible somehow. I tried to change $dh$ at .0005 but this results in an error. Thank you – Api Dec 16 '18 at 20:38
• There is a gap because the gradient is very steep and the steps would need to be smaller as you suggest. The better solution is to use ContourPlot as suggested by @UlrichNeumann. – Bob Hanlon Dec 16 '18 at 21:16
• Hey I was going through this code and I noted that when you change the power in $cdf2[x]$, it does not work anymore. I tried with using as exponent 5 and 10 and got the first time an error and the second time a wrong graph. Do you have any suggestions about how to fix it? Thank you – Api Jan 16 at 16:16

NSolve cannot solve the equation because the equation isn't numeric(depends on h)

For a first insight of the solution use

 ContourPlot[ratio[x, h] == 0.5, {x, -1, 3}, {h, 0, 1},MaxRecursion -> 5, FrameLabel -> {x, h}]


Try

sol[h_?NumericQ] := Values[NSolve[ratio[x, h] == 0.5 && (-5 < x < 5), x]//Flatten]


to get parameter dependent solutions.

sol[.5]
(*{-0.193146, 1.3878, 1.94506}*)
sol[.75]
(*{-0.193147}*)


Unfortunately the solution cannot be plotted, because the number of solutions varies ...

Looking at the Contourplot it is possible to evaluate the Contourline using NDSolve:

First we need one point of the contourline, for example point h==0:

NSolve[{ratio[x, 0] == 1/2, -3 < x < 3}, x, Reals] (*{x -> 2.12759}*)


contourline doesn't change

H = NDSolveValue[{D[ratio[x, h[x]], x] == 0 ,h[2.127591638090098] == 0}, h,{x, -1, 3}]


The left boundary of the solution range

H["Domain"][[1]] (* {-0.193105, 2.23873} *)
x0=%[[1]]


fullfills ratio[x0,h]==1/2 x==x0 is also a solution .

Show[{Plot[H[x], {x, -1, 3}, PlotRange -> {0, 1}],ParametricPlot[{x0, h}, {h, 0, 1}]}, AxesLabel -> {x, h[x]}]


That's it, hope I could help!

• I tried, but then this error was reported: SetDelayed::write: Tag List in {{x->-0.193133},{x->0.774654},{x->2.11263}}[h_?NumericQ] is Protected. – Api Dec 16 '18 at 18:20
• Would you mind to produce the code for plotting it? I am sorry, I am an absolute beginner, I tried some guesses but I always get an empty graph – Api Dec 16 '18 at 18:26
• The code is Plot[Evaluate[ sol[h]], {h, 0, 1}] but evaluation doesn't finish... – Ulrich Neumann Dec 16 '18 at 18:28
• It takes 5 minutes and then it returns a graph which again shows just one of the solutions – Api Dec 16 '18 at 18:34
• Sorry, I don't know why the Plot doesn't work as expected. – Ulrich Neumann Dec 16 '18 at 18:43