# Solving a system of nonlinear equations with parameters

I try to solve a system of two nonlinear first-order conditions for the variables $$\phi_{A}$$ and $$\phi_{B}$$ with two parameters $$\alpha, \beta \in (0,1)$$. Please see example code below. I am looking for symmetric solutions for $$\phi_{A}$$ and $$\phi_{B}$$, which lies between $$0$$ and $$1$$, given it exists.

ClearAll["Global*"]

q1A = (1 - p1A - β + p1B β + α (-1 + p2A + β - p2B β))/((-1 + α^2) (-1 + β^2));
q1B = (1 - p1B - β + p1A β + α (-1 + p2B + β - p2A β))/((-1 + α^2) (-1 + β^2));
q2A = (1 - p2A - β + p2B β + α (-1 + p1A + β - p1B β))/((-1 + α^2) (-1 + β^2));
q2B = (1 - p2B - β + p2A β + α (-1 + p1B + β - p1A β))/((-1 + α^2) (-1 + β^2));

f1 = Simplify[(1 - ϕA)*p1A*q1A + (1 - ϕB)*p1B*q1B];
f2 = Simplify[(1 - ϕA)*p2A*q2A + (1 - ϕB)*p2B*q2B];

{{p1As, p1Bs, p2As, p2Bs}} =
{p1A, p1B, p2A, p2B} /.
Simplify[
Solve[
{D[f1, p1A] == 0, D[f1, p1B] == 0, D[f2, p2A] == 0, D[f2, p2B] == 0},
{p1A, p1B, p2A, p2B}]];

{q1As, q1Bs, q2As, q2Bs} =
Simplify[
{q1A, q1B, q2A, q2B} /. {p1A -> p1As, p1B -> p1Bs, p2A -> p2As, p2B -> p2Bs}];

gA = Simplify[ϕA*(p1As*q1As + p2As*q2As)];
gB = Simplify[ϕB*(p1Bs*q1Bs + p2Bs*q2Bs)];


Problematic command

 Simplify[Solve[{D[gA, ϕA] == 0, D[gB, ϕB] == 0}, {ϕA, ϕB}]]


Simplify[
Solve[{D[gA, ϕA] == 0, D[gB, ϕB] == 0, ϕA == ϕB, ϕA > 0, ϕA < 1}, {ϕA, ϕB}]]


Solution for parameter values $$\alpha=1/2$$ and $$\beta=1/2$$

gA2 = Simplify[gA /. {α -> 1/2, β -> 1/2}];
gB2 = Simplify[gB /. {α -> 1/2, β -> 1/2}];

Simplify[
Solve[{D[gA2, ϕA] == 0, D[gB2, ϕB] == 0, ϕA == ϕB, ϕA > 0, ϕA < 1}, {ϕA, ϕB}]]


Proposed Solution for general α and β

ϕAsol = ((2 - α) (1 - β^2))/(2 - α (1 + β));
Simplify[ϕAsol /. {α -> 1/2, β -> 1/2}]


Solve gives a solution for specific parameter values. For instance, $$\alpha=1/2$$ and $$\beta=1/2$$. However, if I do not specify parameter values for $$\alpha$$ and $$\beta$$ the command does not deliver a solution (I did let it calculate for several hours without result.)

I know that the solution that I am looking for should be $$\phi_A=\phi_B=\frac{\left(2-\alpha\right)\left(1+\beta^{2}\right)}{2-\alpha\left(1+\beta\right)},$$

but unfortunately I cannot receive this result in the above example. Any help is highly appreciated. Thanks in advance. I am new to mathematica and this is my first question on Mathematica StackExchange.

• Welcome to Mathematica.SE, Frank! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. Aug 4, 2019 at 21:07

Substituting your symmetry assumption gives the result you're looking for:

Simplify[Solve[{D[gA, ϕA] == 0, D[gB, ϕB] == 0} /. {ϕA -> ϕ, ϕB -> ϕ}, ϕ]]
(* {{ϕ -> -(((-2 + α) (-1 + β^2))/(-2 + α + α β))}} *)


except the numerator seems to be (2 - α) (1 - β^2) not (2 - α) (1 + β^2).

A numerical example:

α = 0.01;
β = 0.6;
-(((-2 + α) (-1 + β^2))/(-2 + α + α β))
ContourPlot[Evaluate[{D[gA, ϕA] == 0, D[gB, ϕB] == 0}], {ϕA, 0, 1}, {ϕB, 0, 1}]
(* 0.641935 *)
` Looks good.