# Solving parametric function numerically for positive values

So I have a function f1[x,y] that includes a parameter c. I would like to solve the corresponding f1[x,y]==0 for positive x and y (depending on c), with the further restriction that c>=0, but I do not know how to encode all of this properly.

My try was:

Reduce[f[x,y]==0 && x>0 && y>0 &&c>=0, {x,y}]


But it does not seem to work (it is stuck "running"), while playing around a little bit NSolve[f[x,y]==0, {x,y}]/. -> c=1, or any positive c have many solutions, among which some for positive x and y. Conversely, also NSolve[f[x,y]==0 &&x>0 && y>0, {x,y}]/. -> c=1 is stuck running. Am I doing something wrong?

To completely reproduce the framework I am using I would have to add many useless stuff, so I've worked out a minimal reproduction of the problem I have using two functions:

f1[x_, y_] := (50 (2 p^6 (20 + y) - p^2 (-2 + x) x^3 (20 + 3 y) +
p^4 x (3 x (-20 + y) + 2 (60 + y))))/((p^2 + x^2)^3 (20 + y))
f2[x_, y_] := (p^4 (20 + y)^2 + x^4 (20 + y)^2 +
2 p^2 x^2 (1400 - 1000 x + 40 y + y^2))/((p^2 + x^2)^2 (20 + y)^2)

NSolve[f1[x, y] == 0 &&  f2[x, y] == 0, {x, y}] /. p -> 4


Where the NSolve gives as output:

{{x -> 1.01264, y -> -18.804}, {x -> 5.21585,
y -> 24.3412}, {x -> 6.26257, y -> -66.5406}, {x -> 179.176,
y -> -7.27542}, {x -> -4.13331 - 0.549203 I,
y -> -22.4456 + 51.1497 I}, {x -> -4.13331 + 0.549203 I,
y -> -22.4456 - 51.1497 I}, {x -> -0.112599 - 2.93525 I,
y -> 91.3789 + 57.5289 I}, {x -> -0.112599 + 2.93525 I,
y -> 91.3789 - 57.5289 I}, {x -> 0.412582 - 2.45391 I,
y -> -55.0916 - 58.626 I}, {x -> 0.412582 + 2.45391 I,
y -> -55.0916 + 58.626 I}}


Where the second one is clearly positive. However, NSolve[f1[x, y] == 0 && f2[x, y] == 0 && x > 0 && y > 0, {x, y}] /. p -> 4 does not find it and remains stuck. At the same time, Reduce finds roots in the unbounded case and not in the bounded one.

• I should also mention that I do not expect any explicit closed form solutions due to the nature of the function: I'd be happy with positive roots, which do exist (actually, plenty of them according to the first NSolve). – Thanatopseustes Aug 14 '20 at 9:49
• If you reasonably restrict search for solutions with Reduce you should be able to find them. However without definition of f there is nothing to say. If f is too complicated try to provide another example that one might go ahead with. – Artes Aug 14 '20 at 9:50
• Oh, I thought it was a syntax problem (since I am just starting out)... If this is not the case then why NSolve[f[x,y]==0, {x,y}]\.c->1 provides me with many solutions with some positive and NSolve[f[x,y]==0 &&x>0 && y>0, {x,y}]\.c->1 is stuck running? I mean there is no reason for the solutions of the second not to be a subset of the first, if I am not making errors with the coding – Thanatopseustes Aug 14 '20 at 9:58
• Typically you need to present a minimal example so that we can work on it. – yarchik Aug 14 '20 at 11:48
• @Artes Hey hey there's no need to be rude on the internet... I am here precisely because I am a beginner and I wanted some clarification, no need to flex on your superior understanding: I am trying to learn stuff cmon. I've edited the question, hope it sounds "less stupid". – Thanatopseustes Aug 14 '20 at 12:46

NSolve[f1[x, y] == 0 && f2[x, y] == 0 && x > 0 && y > 0, {x, y},
Reals] /. p -> 4


plot:

f1[x_, y_] := (50 (2 p^6 (20 + y) - p^2 (-2 + x) x^3 (20 + 3 y) +
p^4 x (3 x (-20 + y) + 2 (60 + y))))/((p^2 + x^2)^3 (20 + y))
f2[x_, y_] := (p^4 (20 + y)^2 + x^4 (20 + y)^2 +
2 p^2 x^2 (1400 - 1000 x + 40 y + y^2))/((p^2 + x^2)^2 (20 +
y)^2)
ContourPlot[{f1[x, y] == 0, f2[x, y] == 0} /. p -> 4 // Evaluate, {x,
0, 50}, {y, 0, 50}, PlotPoints -> 50]


• It does work like a charm in both the NSolve and Reduce problems... I wonder what happens at the algorithmic level for this to be pivotal. That's the answer though, thanks!. – Thanatopseustes Aug 14 '20 at 13:20

It also works with Solve

Clear["Global*"]

f1[x_, y_] := (50 (2 p^6 (20 + y) - p^2 (-2 + x) x^3 (20 + 3 y) +
p^4 x (3 x (-20 + y) + 2 (60 + y))))/((p^2 + x^2)^3 (20 + y));

f2[x_, y_] := (p^4 (20 + y)^2 + x^4 (20 + y)^2 +
2 p^2 x^2 (1400 - 1000 x + 40 y + y^2))/((p^2 + x^2)^2 (20 + y)^2);

eqns = {f1[x, y] == 0, f2[x, y] == 0, x > 0, y > 0, p >= 0};


The solutions are ConditionalExpressions involving Root expressions

sol[p_] = Solve[eqns, {x, y}];


However, it will evaluate to numeric values for specific values of p meeting the specified conditions.

sol /@ Range[1., 4.]

(* {{{x -> Undefined, y -> Undefined}}, {{x -> 3.05436,
y -> 9.37664}}, {{x -> 4.11556, y -> 17.575}}, {{x -> 5.21585,
y -> 24.3412}}} *)

Plot[Evaluate[{y, x} /. sol[p]], {p, 0, 4},
AxesLabel -> {p, None},
PlotLegends -> Placed[{y, x}, {0.75, 0.4}]] // Quiet
`