# NSolve fails to solve simultaneous equation that has a solution

I'm trying to use Mathematica 12.0 for Students to solve a two equation, two variable simultaneous equation. The context is from microeconomics where 2 firms choose the prices of their respective products p1,p2 to maximise their own profit pi1,pi2 (the equilibrium concept used is Nash). The market shares (q1,q2) is given by the a 'logit' demand system which splits the market based on the prices of the two products.

ClearAll["Global*"];
c1 = 1; c2 = 1; ξ1 = 3; ξ2 = 3; α = 0.1; τ = 0.3;
e1[p1_] := Exp[-α p1 + ξ1];
e2[p2_] := Exp[-α p2 + ξ2]
q1[p1_, p2_] := 100*e1[p1]/(1 + e1[p1] + e2[p2]);
q2[p1_, p2_] := 100*e2[p2]/(1 + e1[p1] + e2[p2]);
pi1[p1_, p2_] := q1[p1, p2] ((1 - τ) p1 - c1);
pi2[p1_, p2_] := q2[p1, p2] ((1 - τ) p2 - c2);
foc1[p1_, p2_] := D[pi1[p1, p2], p1];
foc2[p1_, p2_] := D[pi2[p1, p2], p2]
Simplify[Solve[foc1[p1, p2] == 0, p1]]
Simplify[Solve[foc2[p1, p2] == 0, p2]]
Plot[ArgMax[pi1[p1, p2], p1], {p2, 0, 100}, PlotRange -> Full]
Plot[ArgMax[pi2[p1, p2], p2], {p1, 0, 100}, PlotRange -> Full]
NSolve[{foc1[p1, p2] == 0, foc2[p1, p2] == 0}, {p1, p2}, Reals]


Mathematica gives an expression for Simplify[Solve[foc1[p1, p2] == 0, p1]] (the first order condition for firm 1 choosing p1, taking p2 as given). I did not expect the expression to be clean, so that appears fine.

Solve[foc1[p1, p2] == 0 in this setup should be identical to ArgMax[pi1[p1, p2], so I plot Plot[ArgMax[pi1[p1, p2], p1], {p2, 0, 100}, PlotRange -> Full] to examine the 'best response function' (what firm1 should set p1 to given a range of values for p2). The intersection of the 2 firms' best response functions gives the equilibrium. The best response functions appears to have good properties that should allow its intersections to be obtain relatively easily. In this example, the firms are symmetric, so the solution to the intersection of is where the 45 degree line cuts the plots, so it does have a solution (I believe the setup should mean it is unique).

However, NSolve[{foc1[p1, p2] == 0, foc2[p1, p2] == 0}, {p1, p2}, Reals] does not output an answer.

Have I done something wrong with my setup or is there anything I can specify for Mathematica to get it to solve the simultaneous equation?

• Since you are OK with a numerical solution, and you should have a good starting point values from your plots, try using FindRoot instead of NSolve: for instance, FindRoot[{foc1[p1, p2] == 0, foc2[p1, p2] == 0}, {{p1, 5}, {p2, 5}}] returns a solution very quickly. Dec 23 '19 at 17:35
• In Documentation is written: "NSolve deals primarily with linear and polynomial equations". Dec 23 '19 at 17:55
• Plot the functions to get an impresson where they both are zero. Plot3D[{0, foc1[p1, p2], foc2[p1, p2]}, {p1, 0, 1000}, {p2, 0, 1000}, PlotStyle -> {Yellow, Blue, Red}] . Dec 23 '19 at 19:19

The problem can be solved by exploiting the symmetry of the equations.

Clear["Global*"]


Use exact values for the constants to enable an exact solution.

EDIT: corrected typo in value for α (Thanks to Akku14)

c1 = 1; c2 = 1; ξ1 = 3; ξ2 = 3; α = 1/10; τ = 3/10;
e1[p1_] := Exp[-α p1 + ξ1];
e2[p2_] := Exp[-α p2 + ξ2]
q1[p1_, p2_] := 100*e1[p1]/(1 + e1[p1] + e2[p2]);
q2[p1_, p2_] := 100*e2[p2]/(1 + e1[p1] + e2[p2]);
pi1[p1_, p2_] := q1[p1, p2] ((1 - τ) p1 - c1);
pi2[p1_, p2_] := q2[p1, p2] ((1 - τ) p2 - c2);
foc1[p1_, p2_] := D[pi1[p1, p2], p1];
foc2[p1_, p2_] := D[pi2[p1, p2], p2];

Plot3D[Evaluate@{0, foc1[p1, p2], foc2[p1, p2]},
{p1, 0, 50}, {p2, 0, 50},
ClippingStyle -> None]


Simplify the equations assuming that {p1 > 0, p2 > 0}

eqns = foc1[p1, p2] == 0 && foc2[p1, p2] == 0 //
Simplify[#, {p1 > 0, p2 > 0}] &

(* E^(p1/10) (E^3 + E^(p2/10)) (-80 + 7 p1) == 70 E^(3 + p2/10) &&
E^(p2/10) (E^3 + E^(p1/10)) (-80 + 7 p2) == 70 E^(3 + p1/10) *)


From the symmetry of the plot and the equations, assume that p1 == p2

The exact solution is a Root object

(sol = p1 /.
Solve[{eqns /. p2 -> p1 // Simplify, p1 > 0}, p1][[1]]) // InputForm

(* Root[{-150 + 7*#1 + E^(#1/10)*
(-80/E^3 + (7*#1)/E^3) & , 18.\
94202035330226284354267930111274876\
86220.602059668174107}] *)


As indicated in the Root expression, this is approximately,

sol // N[#, 20] &

(* 18.942020353302262844 *)


Verifying that this value for both p1 and p2 satisfies the equations

eqns /. N[{p1 -> sol, p2 -> sol}, 20]

(* True *)


Similarly for NSolve

solN = p1 /. NSolve[
{eqns /. p2 -> p1 // Simplify, p1 > 0}, p1,
WorkingPrecision -> 20][[1]]

(* 18.942020353302262844 *)

eqns /. {p1 -> solN, p2 -> solN}

(* True *)

• A little typo? The OP askes for alpha = 1/10 . Therefore sol // N[#, 20] & is (* 18.942020353302262844 *)` . Dec 23 '19 at 22:50
• @Akku14 - Thanks. Corrected. Dec 23 '19 at 23:08