Solving a system of nonlinear equations with parameters

Question

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}]]

Alternative command with additional restrictions

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.

Answer 1

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.