How to solve this system of multivariate polynomial equations for $0<x_7<x_6<x_8 \le 1$? Groebner basis maybe?

Question

I am reformulating my question according to the guidelines I was given. I have the following problem: I cannot find a way to solve the system of equations further down. This is the calculations from top to bottom. All values are in $\mathbb{R}$. \begin{align*} a_1 &:=\frac{x_3-x_1}{2x_7}+\frac{x_7}{2}\\ a_2 &:=\frac{x_2-x_3}{2(x_6-x_7)}+\frac{x_7+x_6}{2}\\ a_3 &:=\frac{x_4-x_2}{2(x_8-x_6)}+\frac{x_8+x_6}{2} \end{align*} We have: \begin{align*} \pi_A &= x_1\cdot a_1 + x_2\cdot (a_3-a_2)\\ \pi_B &= x_3\cdot (a_2-a_1) + x_4\cdot (1-a_3) \end{align*} I then maximize $\pi_A$ w.r.t $x_1$ and $x_2$ and $\pi_B$ w.r.t $x_3$ and $x_4$. So: \begin{align*} \frac{\partial \pi_A}{\partial x_i}&\stackrel{!}{=} 0 \quad i \in \{1,2\}\\ \frac{\partial \pi_B}{\partial x_i}&\stackrel{!}{=} 0 \quad i \in \{3,4\} \end{align*} We solve for $x_i$ for $i= 1\ldots 4$. The solutions are unique and we substitute them back into our functions $\pi_A$ and $\pi_B$. Now this is where my problems and my questions start. I now want to solve: \begin{align} \frac{\partial \pi_A}{\partial x_6}& = 0 \\ \frac{\partial \pi_B}{\partial x_7}& = 0 \\ \frac{\partial \pi_B}{\partial x_8}& = 0 \end{align} I am only interested in solutions that satisfy $0<x_7<x_6<x_8 \le 1$. Any ideas maybe? I am a little familiar with Groebner bases and Buchbgerger's algorithm but have not succeeded trying to apply them. Also here a link to files in Maple and Matlab that I have created in trying to solve this: https://github.com/fabsongithub/Interlacing

For the people interested in the background of those calculations. I use them to find conclusions about a spatial competition model on the Hotelling line.

Answer 1

Using Mathematica, once we have the desired set of equations:

Clear[x1, x2, x3, x4];

a1 = (x3 - x1)/(2 x7) + x7/2;
a2 = (x2 - x3)/(2 (x6 - x7)) + (x7 + x6)/2;
a3 = (x4 - x2)/(2 (x8 - x6)) + (x8 + x6)/2;

πA = x1 a1 + x2 (a3 - a2);
πB = x3 (a2 - a1) + x4 (1 - a3);

sol = Solve[{D[πA, x1] == 0, D[πA, x2] == 0,
             D[πB, x3] == 0, D[πB, x4] == 0},
            {x1, x2, x3, x4}];
{{x1, x2, x3, x4}} = FullSimplify[sol[[All, All, 2]]];

eqns = FullSimplify[{D[πA, x6], D[πB, x7], D[πB, x8]}
        18 ((x6 + x7)^2 - 4 x6 x8)^3/{x7 - x8, x6, 1}];

we can obtain a numerical solution as follows:

NSolve[{eqns == {0, 0, 0}, 0 < x7 < x6 < x8 <= 1}]

{}

which means that there are no acceptable solutions, as can also be ascertained graphically:

ContourPlot3D[{eqns[[1]] == 0, eqns[[2]] == 0, eqns[[3]] == 0},
              {x6, 0, 1}, {x7, 0, 1}, {x8, 0, 1}, AxesLabel -> {"x6", "x7", "x8"},
              PlotPoints -> 50, RegionBoundaryStyle -> None,
              RegionFunction -> Function[{x6, x7, x8}, 0 < x7 < x6 < x8 <= 1]]

$\quad\quad\quad\quad\quad\quad$

where of course this is appreciable only in Mathematica by rotating the 3D graph.