Why is this equation difficult to solve numerically?

I'm solving a fairly simple equation :

w[p1_, p2_, xT_] :=
94.8*cv*p1*y[(p1 - p2)/p1, xT]*Sqrt[(p1 - p2)/p1*mw/t1];

y[r_, xT_] := 1 - (1.4 r)/(3 xT*γ) /. γ -> 1.28;

sol = NSolve[{w[p1, p2, 0.66] == 30, p2 == 1.07}, {p1, p2}] /. {cv ->
1.77, t1 -> 318, mw -> 38};


Mathematica has no problem solving this case. However, when I change the y function to:

y[r_, xT_] :=
Max[0.667, 1 - (1.4 r)/(3 xT*γ)] /. γ -> 1.28


Then the calculation does not seem to converge. Is there a reason why this should be difficult ? Can it be implemented differently ?

• If probably fails because Max plus the nonnumeric parameters are beyond the methods of NSolve, which has a limited scope (see docs). Works if you substitute the values of the parameters first: NSolve[{w[p1, p2, 0.66] == 30, p2 == 1.07} /. {cv -> 1.77, t1 -> 318, mw -> 38}, {p1, p2}] – Michael E2 Aug 24 '16 at 10:21
• OK, so in general it's a good idea to replace the coefficients with their numerical values as early as possible. – Whelp Aug 24 '16 at 11:42
• Yes, for NSolve, as I believe some methods of NSolve work only (or better?) on determinate problems. – Michael E2 Aug 24 '16 at 12:16

The discontinuity makes it so that

NSolve was unable to solve the system with inexact coefficients

The solution of NSolve[] in this case is to:

The answer was obtained by solving a corresponding exact system and numericizing the result.

However, it is impossible to do so with your strategy of list replacements.

cv = 1.77; t1 = 318; mw = 38; γ = 1.28;
w[p1_, p2_, xT_] :=
94.8*cv*p1*y[(p1 - p2)/p1, xT]*Sqrt[(p1 - p2)/p1*mw/t1];
y[r_, xT_] :=
Max[0.667, 1 - (1.4 r)/(3 xT*γ) /. γ -> 1.28];
sol = NSolve[{w[p1, p2, 0.66] == 30, p2 == 1.07}, {p1, p2}]


{{p1 -> 1.32278, p2 -> 1.07}}

Which corresponds to the third answer obtained by the original problem. Note that the numericization of an exact system doesn't leave imaginary rounding errors like with the original problem, don't forget to Chop[].

My original answer was a misinterpretation. I post this in case it is helpful.

y[r_, xT_, g_] := 1 - (1.4 r)/(3 xT*g)
w[p1_, p2_, xT_, cv_, mw_, t1_] :=
94.8*cv*p1*y[(p1 - p2)/p1, xT, 1.28]*Sqrt[(p1 - p2)/p1*mw/t1];
ans = p2 /. NSolve[w[1.07, p2, 0.66, 1.77, 38, 318] == 30, p2]
plot1 = Plot[{w[1.07, p2, 0.66, 1.77, 38, 318],
10 y[(1.07 - p2)/1.07, 0.66, 1.28]}, {p2, 0, 1},
MeshFunctions -> (#2 &), Mesh -> {{30}},
MeshStyle -> {Red, PointSize[0.02]}, GridLines -> {ans, {30}},
Frame -> True,
FrameTicks -> {{Automatic,
Table[{j, j/10.}, {j, 0, 40, 5}]}, {Automatic, None}},
FrameLabel -> {"p1", "w", "", "y"}, ImageSize -> 300,
PlotRange -> {0, 40}]
y2[r_, xT_, g_] := Max[0.667, 1 - (1.4 r)/(3 xT*g)]
w2[p1_, p2_, xT_, cv_, mw_, t1_] :=
94.8*cv*p1*y2[(p1 - p2)/p1, xT, 1.28]*Sqrt[(p1 - p2)/p1*mw/t1];
ans2 = p2 /. NSolve[w2[1.07, p2, 0.66, 1.77, 38, 318] == 30, p2]
plot2 = Plot[{w2[1.07, p2, 0.66, 1.77, 38, 318],
10 y2[(1.07 - p2)/1.07, 0.66, 1.28]}, {p2, 0, 1},
MeshFunctions -> (#2 &), Mesh -> {{30}},
MeshStyle -> {Red, PointSize[0.02]}, GridLines -> {ans2, {30}},
Frame -> True,
FrameTicks -> {{Automatic,
Table[{j, j/10.}, {j, 0, 40, 5}]}, {Automatic, None}},
FrameLabel -> {"p1", "w", "", "y"}, PlotRange -> {0, 40},
ImageSize -> 300]


Visualizing:

• just avoided all the substitions
• always find plotting useful
• desired solution is selected (first form->{0.676056, 0.134129}, second case->0.676056)