The "short answer" is as follows.
The point $P$ is in general not constructible, for instance, for the triangle with sides $6$, $9$, $13$
constructing $P$ implies we can solve a polynomial equation $\Pi(K)=0$, where
$\Pi=a_7K^7+a_6K^6+a_5K^5+a_4K^4+a_3K^3+a_2K^2 + a_1K+a_0$
is an irreducible polynomial of degree seven with rational coefficients,
and a root $K$ of it is the square of the common distance of $P$ to the vertices.
I will use $K$ for both the indeterminate when writing $\Pi$,
and also the particular root value(s) $K\in\overline{\Bbb Q}$ with corresponding to real positive numbers.
A simpler (particular) case is when $\Delta ABC$ is isosceles.
Then $\Pi$ factorizes as $\Pi=\Pi_2^2\cdot\Pi_3$ with $\deg \Pi_2=2$, $\deg \Pi_3=3$,
$\Pi_2$ without roots in $\Bbb R$.
So the possible $K$ values come from the real positive roots of $\Pi_3$. In explicit cases we can easier check $\Pi_3$
has no rational roots. So the constructibility fails.
The point $P$ fails also to be unique in the plane - as stated by the OP - in infinitely many cases.
To see which cases lead to a failure,
explicit computations of the coefficients of $\Pi$ in terms of the parameters $a,b,c$ are needed.
It turns out that in "most cases" there is only one sign change in the list of coefficients of $\Pi$,
the first coefficients $a_7$, $\color{blue}{a_6}$, $a_5$, $a_4$, $a_3$ being positive, and $a_2$, $a_1$, $a_0$ being negative,
so Decartes' Rule of Signs insures an unique $K$.
However, in "few cases" the coefficient $\color{blue}{a_6}=\color{blue}{-p^4} + 12 p^2 S^2 E + 12 S^4 E^2 + 192 S^6$
is negative, Decartes' Rule of Signs predicts at most three positive roots, and indeed, we have three such roots $K>0$,
leading to three solutions.
Here $p=abc$, $E=a^2+b^2+c^2$, and $S$ is the area of $\Delta ABC$.
Counterexamples should thus have a "small area" $S$, but "big product" $p$.
It is natural to search then for a counterexample with $c$ being "almost" $a+b$, and this produces them quickly, as checked by computer.
We have however a uniqueness of $P$ in the interior of the given triangle.
This, together with the existence, may be seen by geometric arguments as in the other answers.
For the existence of a solution i have a deformation argument.
For existence and uniqueness in the interior this answer gives some hints / details for the deformation.
(Making things concrete goes beyond the question of the OP.)
Given this verdict, a / my proof cannot work by geometric means,
so analytic tools are needed, my choice was to use barycentric coordinates.
I will try to follow in presentation a minimal path, however providing the needed details
makes minimal not really compact.
The question raised a big echo, it is indeed an interesting question.
(I'm afraid this answer shows things are not so simple, but please give it a chance.)
To control the computations, i need to use computer algebra support (CAS), my choice of weapons is sage,
it also has a great merit in making the code pretty readable for a mathematician, since the methods used
are named to be easily digested by a mathematician.
As reference for barycentric coordinates i will use
Barycentric Coordinates for the Impatient, Max Schindler, Evan Chen, July 13, 2012
The detailed answer starts now.
Let $a,b,c$ be the lengths of the sides of a general non-degenerated triangle $\Delta ABC$.
Let $x,y,z$ with $1=x+y+z$ be the coordinates of the unknown point $P$. So $P=(x,y,z)$ in notation.
Then $A,B,C;P$ and the points $D=AP\cap BC$, $E=BP\cap CA$, $F=CP\cap AB$.
have explicit barycentric descriptions:
$$
\begin{aligned}
A &= (1,0,0)\ ,\\
B &= (0,1,0)\ ,\\
C &= (0,0,1)\ ,\\[2mm]
P &= (x,y,z)\ ,\\
D &= [0:y:z]=\left(0,\frac {y}{y+z},\frac {z}{y+z}\right)\ ,\\
E &= [x:0:z]=\left(\frac {x}{x+z},0,\frac {z}{x+z}\right)\ ,\\
F &= [x:y:0]=\left(\frac {x}{x+y},\frac {y}{x+y},0\right)\ ,\\[3mm]
&\qquad\text{and corresponding displacement vectors are}\\[3mm]
\overrightarrow{DP} &= D - P = \left(-x,\frac {xy}{y+z},\frac {xz}{y+z}\right) = \frac x{y+z}(-(y+z),y,z) \ ,\\
\overrightarrow{EP} &= E - P = \left(\frac {xy}{x+z},-y,\frac {yz}{x+z}\right) = \frac y{z+x}(x,-(z+x),z) \ ,\\
\overrightarrow{FP} &= F - P = \left(\frac {xz}{x+y},\frac {yz}{x+y},-z\right) = \frac z{x+y}(x,y,-(x+y)) \ ,\\[3mm]
&\qquad\text{and corresponding squared lengths are}\\[3mm]
|DP|^2 &= \frac{x^2}{(y+z)^2}\Big(\ -a^2yz + (b^2z+c^2y)(y+z)\ \Big) = \frac{x^2}{(1-x)^2}\Big(\ Q + b^2z + c^2y\ \Big)\ ,\\
|EP|^2 &= \frac{y^2}{(z+x)^2}\Big(\ -b^2zx + (c^2x+a^2z)(z+x)\ \Big) = \frac{y^2}{(1-y)^2}\Big(\ Q + c^2x + a^2z\ \Big)\ ,\\
|FP|^2 &= \frac{z^2}{(x+y)^2}\Big(\ -c^2xy + (a^2y+b^2x)(x+y)\ \Big) = \frac{z^2}{(1-z)^2}\Big(\ Q + a^2y + b^2x\ \Big)\ ,\\
&\qquad\text{where}\\
Q &= -a^2 yz -b^2 zx -c^2xy\ .
\end{aligned}
$$
Let $K>0$ be the common value of the squared distances from $P$
to each of the points $D,E,F$, i.e. $K=|DP|^2=|EP|^2=|FP|^2$.
Then we have to solve the following system of equations in the unknowns $x,y,z;K,Q\in \Bbb R$, $K>0$:
$$
\tag{$\dagger$}
$$
$$
\left\{
\begin{aligned}
1 &= x+y+z\ ,\\
Q &= -a^2yz -b^2zx -c^2xy\ ,\\[2mm]
K(1-x)^2 &= x^2(Q + b^2 z + c^2 y)\ ,\\
K(1-y)^2 &= y^2(Q + c^2 x + a^2 z)\ ,\\
K(1-z)^2 &= z^2(Q + a^2 y + b^2 x)\ ,
\end{aligned}
\right.
$$
(Solutions were introduced by eliminating denomiators, for instance $A=(1,0,0)$ is now a solution.)
Let $J$ be the ideal generated by the above equations (rewritten in terms vanishing of expressions).
From now on, we will work algebraically with the system $(\dagger)$ above.
(Since the parameters $a,b,c$ appear only through their squares, it may be convenient to
use short hand notations $A=a^2$, $B=b^2$, $C=c^2$ for them.)
(There should be no confusion with the vertices
$A,B,C$ of the given triangle.)
Now i am trying to address the points from the question in terms of this system.
Constructibility of $P$, special case of the $6$, $9$, $13$ triangle:
This triangle is used for the purpose of (numerically) searching for knwon (general) centers in a triangle, see also
ETC, Search_6_9_13 .
If $P$ is somehow constructible (starting from some fixed values $a,b,c$), then $K$ is also constructible (by rule and compass constructions),
so starting with $6,9,13\in\Bbb Q$
one should obtain that $K$ is an algebraic number in an extension of degree among $1,2,4,8,16,\dots$ -
however, it turns out that $K$ is the root of an irreducible polynomial $\Pi$ of degree seven. The sage code computing $\Pi$
is postponed.
In this special case:
$$
\tag{$1$}
\Pi =
K^7
+ \frac{2420849677}{40884480}\; K^6
+ \frac{2635885}{2106 }\; K^5
+ \frac{920646335}{82134 }\; K^4
+ \frac{349438600}{9477 }\; K^3
- \frac{10186414000}{123201 }\; K^2
- \frac{8451520000}{9477 }\; K
- \frac{215129600000}{123201}
\ .
$$
So the constructibility fails in this case, so it fails.
Let us observe that there is exactly one sign change in the coefficients of $\Pi$.
This happens "often", as seen in the next section.
Below there will be an other example ($61,61,120$ isosceles triangle) where the
above polynomial splits as a product of a squared quadratic and a cubic polynomial,
the cubic part only has real roots, we have an explicit formula for them, and it is easy to check
there is no rational root.
The minimal polynomial of $K$ in the general case.
In the code part, a polynomial $\Pi$ of degree seven in $R[K]$ over the ring $R=\Bbb Q[a,b,c]$ is computed, so that
the special value $K=|DP|^2=|EP|^2=|FP|^2$ (for each solution $P$ of our problem) satisfies $\Pi(K)=0$.
Explicitly,
$$
\tag{$2$}
$$
$$
\begin{aligned}
\Pi &=
a_7 K^7 +
a_6 K^6 +
a_5 K^5 +
a_4 K^4 +
a_3 K^3 +
a_2 K^2 +
a_1 K +
a_0\ ,\\[3mm]
&\qquad\text{where}\\[3mm]
a_7 &= 64 S^2 p^2\ ,\\
a_6 &= \color{blue}{-p^4} + 12 p^2 S^2 E + 12 S^4 E^2 + 192 S^6\ ,\\
a_5 &= 2 S^4 \;(E^3 + 13p^2 + 34 S^2 E)\ ,\\
a_4 &= S^4 \;(p^2E + 10 S^2 E^2 + 103 S^4)\ ,\\
a_3 &= S^6 \;(2p^2 + 17 S^2 E)\ ,\\
a_2 &= -\frac 14 S^8 \;(E^2 - 40 S^2)\ ,\\
a_1 &= -E \;S^{10} \ ,\\
a_0 &= -S^{12}\ .
\end{aligned}
$$
The coefficients are homogeneous of degrees $\deg a_7=10$, $\deg a_6=12$, ... , $\deg a_0=24$,
when considering the weights $\deg S=\deg E=\deg K=2$, $\deg p=3$.
If $K$ is considered with $\deg K=2$, then $\Pi$ is homogeneous of degree $24$.
These explicit formulas allow to control the signs of the coefficients of $\Pi=\Pi(K)$.
We have (at least) one sign change obtained when passing from $a_3$ to $a_2$.
To see this note that
$$
E^2 - 40 S^2 > E^2 - 48 S^2 = 4(a^4+b^4+c^4-a^2b^2-b^2c^2-c^2a^2)>0\ .
$$
An other change of signs may occur only around $a_6$.
In case $a_6>0$ there is exactly one change of signs at all,
so by Decartes' rule of signs we have exactly one real positive root.
So uniqueness of $K$ is insured in case of $a_6>0$.
Let us analyze this in detail.
Further relations:
Let $x,y,z;K$ be the solutions of the system $(\dagger)$. Then
there is a relation joining $K$ with each of the variables $x,y,z$:
$$
\tag{$3$}
$$
$$
\begin{aligned}
0 &= 4S^2 x^3 - 2S^2x^2 +K(b^2+c^2-a^2)x^2 + 8 xK^2 - 2K^2 \ ,\\
0 &= 4S^2 y^3 - 2S^2y^2 +K(c^2+a^2-b^2)y^2 + 8 yK^2 - 2K^2 \ ,\\
0 &= 4S^2 z^3 - 2S^2z^2 +K(a^2+b^2-c^2)z^2 + 8 zK^2 - 2K^2 \ ,
\end{aligned}
$$
and relations for $(y,K)$ and $(z,K)$ can be written analogously.
We also have:
$$
\tag{$4$}
$$
$$
\begin{aligned}
0 &= a^2\;y^2z^2 + 4K\;yz(y+z) - K(y+z)^2\ ,\\
0 &= b^2\;z^2x^2 + 4K\;zx(z+x) - K(z+x)^2\ ,\\
0 &= c^2\;x^2y^2 + 4K\;xy(x+y) - K(x+y)^2\ ,
\end{aligned}
$$
as found in the code section #4.
(These relations can be used maybe in some way to show algebraically the existence
and uniqueness in the interior of $\Delta ABC$
for the point matching a fixed $K$, solution of the above polynomial $\Pi$.
I tried to make it, but it made me tired.)
Uniqueness of $P$ fails in the plane:
We have uniqueness,
in case there is a unique positive root $K$
of $\Pi$.
For instance, in the "often" case of a coefficient $a_6>0$ of $\Pi$,
we know from above that there exists exactly one $K>0$ so that for a $P$ with $|PD|^2=|PE|^2=|PF|^2$ (in case of its existence)
the common value of the above squared distances is $K$.
However cases can be found where such a $K$ is not unique, and we can check that multiple solutions arise.
In the code section #5 there is the following solution found for the triangle with sides $a=b=61$, $c=120$.
The height $h_c$ is $11$. The solutions are roughly:
$$
\begin{aligned}
x = y & = -1.749464267458191 \dots & z &= 4.498928534916382 \dots & K &= 2449.08331\dots\\
x = y & = -0.5905364374782210\dots & z &= 2.181072874956442 \dots & K &= 575.606545\dots\\
x = y & = 0.1209924404736021\dots & z &= 0.7580151190527959\dots & K &= 69.5250174\dots\\
\end{aligned}
$$
Only the last point lives inside $\Delta ABC$.
$K$ is a root of the polynomial $(121K^3 - 374400K^2 + 196020000 K - 11859210000)$.
(For the point $P$ inside $\Delta ABC$ the corresponding $F$ is the mid point of $AB$,
the height $CF$ has length $11$ and $PF$ is approximatively $\sqrt{69.5250174\dots}$. The other two points
have the same $F$. Using this point as origin,
$$\overrightarrow{FP}=x(\overrightarrow{FA}+\overrightarrow{FB})+z\overrightarrow{FC}
=z\overrightarrow{FC}
$$
in all three listed cases. And indeed, $121\;z^2=|CF|^2 \; z^2=|PF|^2=K$ matches the values for $K$
for each $z$ in the list.
For these $x,y,z$ we can also write down exact polynomials having them as roots.
The three $x$-values are for instance roots of the polynomial
$\displaystyle x^3 + \frac{537}{242} x^2 + \frac 34 x - \frac 18$.
none of them is rational.
Of course, "small deformations" lead also to situations with three solutions.
For the existence, geometric ideas may work better.
Existence and uniqueness of $P$ in the interior of $\Delta ABC$:
So we restrict the values for the barycentric coordinates $x,y,z$ to positive values,
$x, y, z>0$, $x+y+z=1$.
Consider the following picture:
(We do not know right now that an interior $P^*$ exists so that for the cevians $AD^*$, $BE^*$, $CF^*$ through $P^*$ we have $|P^*D^*|=|P^*E^*|=|P^*F^*|$.
Also its uniqueness is still an issue.)
(Also the other answers suggest a possibility to "get closer to the / a solution $P^*$ - if it exists.)
I propose to use two geometrical schemes to "get" (a) $P^*$ in the limit.
First scheme:
Start with a point $P=P_0$. Compare the distances from $P$ to $D,E,F$ and take the maximal one, say it is $PE$.
Let a point $Q=P_1$ slide from $P$ to $E$.
Then we stop when
$f_P(Q):=2|QE_1|-|QD_1|-|QF_1|$
vanishes first $f_P(Q)=0$.
Then iterate.
For the convergence we need some argument. (Is this procedure / function $P\to Q$ a contraction?)
Second scheme:
Start with an interior point $P=(x,y,z)$, and consider the squared distances $PD^2$, $PE^2$, $PF^2$.
Then consider the point $Q=[X:Y:Z]$ where $X,Y,Z$ are weighted versions of $(x,y,z)$ as follows:
$$
\begin{aligned}
X &=x\cdot \frac{PE^2+PF^2}{2(PD^2+PE^2+PF^2)}\ ,\\
Y &=y\cdot \frac{PF^2+PD^2}{2(PD^2+PE^2+PF^2)}\ ,\\
Z &=z\cdot \frac{PD^2+PE^2}{2(PD^2+PE^2+PF^2)}\ .
\end{aligned}
$$
(Intuition: If for instance $PE^2$ is maximal among $PD^2$, $PE^2$, $PF^2$, then we would like to move $P$ along $PE$ towards $E$.
This corresponds to making $y$ smaller. We can try to multiply $y$ with some "controlled" subunitary factor - and
$\frac{PF^2+PD^2}{2(PD^2+PE^2+PF^2)}$ is such a factor.
At the end, the point $[X:Y:Z]$ has to be normed again, so the result is $Q=\frac 1{X+Y+Z}(X,Y,Z)$.
Of course, this is only an Ansatz, one has to show it works.
I only have numerical support so far - see the code section #6.)
Third scheme:
This is less explicit, but we do not have issues as the ones described above.
Let $K\ge0$ be a parameter, we let it variate from $0$ to $\infty$. For each such $K$ we draw three curves.
The $D$-curve is the locus of all points $P$ obtained as follows. Let $D$ slide along the line $BC$.
Draw the cevian $AD$. Let $P$ be on the ray $[DA$ such that $PD^2=K$. Then the $D$-curve is the locus of the points $P$ as $D$ runs on $BC$.
The $D$-curve has the line $BC$ as double asymptote, it is symmetric w.r.t. the height from $A$ and for $K\to 0_+$ it tends to $BC$.
Use as orientation for it the direction that in the limit corresponds to the direction of $BC$ from $B$ to $C$.
Similarly consider the $E$-curve, and the $F$-curve.
They are oriented corresponding to the directions that in the limit give the direction of $CA$ from $C$ to $A$, and of $AB$ from $A$ to $B$.
For $K\to 0_+$, in the limit $K=0$, the three curves become $BC$, $CA$, $AB$, and the oriented area between these curves
is $S>0$. Now let $K$ grow. At some point, e.g. when $K$ is greater than the square of the biggest height, the oriented area becomes negative.
So at some point $K^*$ this area is zero. This corresponds to the the case when the three curves are passing through one point.
(I do not have a simple argument to show this is an interior point.)
This elucidates the existence.
$\square$
A rough picture for the three curves at their intersection in the zoom of our objective is as follows:
I hope it is clear how the three "lens" of the zoom are moving. To fix idea, assume $\Delta ABC$ has all angles $<90^\circ$,
so the orthocenter $H$ is in its interior. Let $AA'$, $BB'$, $CC'$ be the heights, intersecting in $H$.
- For $K=0$ the area between the three curves is $\Delta ABC$.
- Arrange that $HA'\le HB'\le HC'$ by permuting $A,B,C$.
- when $K=A'H^2$, we still have a positive area between the three curves.
- when $K=C'H^2$, we have a negative area between the three curves.
- so the value $K^*$ is squared between these values.
In case of an obtuse angle - say in $A$ - use an alternative zooming objective
with mobile lens for the $E$- and $F$-curves and keep the third lens constant to be the line $BC$ (instead of the $D$-curve)
to see that the intersection of the three curves still has to be an interior point.
Computer algebra support.
Code #1: The value of $K=|PD|^2 =|PD|^2 =|PD|^2$ in the case of the $6,9,13$ triangle.
a, b, c = 6, 9, 13
R.<x,y,z,K> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g}\n')
And we obtain:
Generator(s) of the elimination ideal JK after eliminating x, y, z:
K^9 + 2420849677/40884480*K^8 + 2635885/2106*K^7 + 920646335/82134*K^6 + 349438600/9477*K^5
- 10186414000/123201*K^4 - 8451520000/9477*K^3 - 215129600000/123201*K^2
It turns out that the above polynomial is $K^2$ times some irreducible polynomial of degree seven.
To have the entry that should be compared with ETC, we ask for the points in the ring of real algebraic numbers:
sage: J.variety(ring=AA)
[{K: 0, z: 0, y: 0, x: 1},
{K: 0, z: 0, y: 1, x: 0},
{K: 0, z: 1, y: 0, x: 0},
{K: 3.973344192056688?,
z: 0.5361103522937736?, y: 0.2667643469973961?, x: 0.1971253007088303?}]
(New roots were introduced after multiplication with denominators.)
Only the last real point is significant. In ETC we have the value for
$2\cdot\operatorname{Area}(\Delta ABC)\cdot \frac xa$. The area of the $6,9,13$ triangle is $\sqrt{14(14-6)(14-9)(14-13)}=\sqrt{560}$.
So we try to match the following value with the 6-9-13-search:
sage: x0 = J.variety(ring=AA)[-1][x]
sage: etc_match = 2 * sqrt(560) * x0/6
sage: etc_match.n(200)
1.5549453416812577768807502448833414277833445395269911652987
This was already done in the comment of
Peter Taylor. To have again a convincing statement against constructibility,
the value x0 used above is an algebraic number, root of an irreducible polynomial in $\Bbb Q[x]$
of degree seven:
sage: x0.minpoly()
x^7 - 28472791/10221120*x^6 + 46797133/15331680*x^5 - 52400213/30663360*x^4
+ 28733/54756*x^3 - 5173/54756*x^2 + 280/13689*x - 35/13689
Code #2: The value of $K=|PD|^2 =|PD|^2 =|PD|^2$ in the case of an isosceles, general triangle with sides $a,a,b$.
To avoid the factor $K^2$, i will assume $K$ invertible below.
R.<a,b, x,y,z, K,K_inv> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, a, b, x, y, z, K),
eq(a, b, a, y, z, x, K),
eq(b, a, a, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z, K_inv])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g.factor()}\n')
And we obtain:
Generator(s) of the elimination ideal JK after eliminating x, y, z:
(1/256) * (16*a^4*b^2 - 8*a^2*b^4 + b^6 + 64*a^4*K + 32*a^2*b^2*K - 12*b^4*K + 256*a^2*K^2)
* (16*a^4*b^4 - 8*a^2*b^6 + b^8 - 32*a^2*b^4*K + 8*b^6*K
- 1024*a^2*b^2*K^2 + 272*b^4*K^2 - 4096*a^2*K^3 + 1024*b^2*K^3)
(Result was manually rearranged.)
Even in this case, the constructibility fails. For instance for the sides $2,2,3$ we have:
a, b, c = 2, 2, 3
R.<x,y,z, K,K_inv> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z, K_inv])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g.factor()}\n')
This gives:
Generator(s) of the elimination ideal JK after eliminating x, y, z:
(1/7340032) * (1024*K^2 + 1204*K + 441) * (7168*K^3 + 14832*K^2 + 4536*K - 3969)
The quadratic factor has no real roots.
The cubic factor has exactly one positive root,
$$
\frac3{448}
\left( -103 + \sqrt[3]{5(111205 + 6272\sqrt{105}} + \sqrt[3]{5(111205 - 6272\sqrt{105}} \right)
\\ \approx 0.3643855138832992\dots
$$
Constructibility fails.
Code #3:
We consider $K=|PD|^2 =|PE|^2 =|PF|^2$ (if the needed $P$ exists) as a function of the sides $a,b,c$
of a triangle. Then there is a polynomial $\Pi\in\Bbb Q[a,b,c]\;[\kappa]$ of degree seven w.r.t. $\kappa$, which vanishes in $K$, $\Pi(a,b,c;K)=0$.
Its coefficients, expressed in terms of the quantities $E = a^2+b^2+c^2$, $S^2=s(s-a)(s-b)(s-c)$ (squared area), and $p=abc$ (product) are obtained as follows:
R.<a,b,c, x,y,z, K,K_inverse, p,SS,E> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1,
K * K_inverse - 1,
SS - 1/16 * (a+b+c) * (a+b-c) * (b+c-a) * (c+a-b),
E - (a^2 + b^2 + c^2),
p - a*b*c,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
for g in J.elimination_ideal([x, y, z, K_inverse, a, b, c]).groebner_basis():
for k in range(g.degree(K), 0, -1):
print(f'Coefficient of K^{k} is {g.coefficient(K**k).factor()}')
print(f'Coefficient of K^0 is {g.subs({K : 0}).factor()}')
And the coefficient are explicitly:
Coefficient of K^7 is (-64) * SS * p^2
Coefficient of K^6 is (-1) * (-p^4 + 12*p^2*SS*E + 12*SS^2*E^2 + 192*SS^3)
Coefficient of K^5 is (-2) * SS^2 * (E^3 + 13*p^2 + 34*SS*E)
Coefficient of K^4 is (-1) * SS^2 * (p^2*E + 10*SS*E^2 + 103*SS^2)
Coefficient of K^3 is (-1) * SS^3 * (2*p^2 + 17*SS*E)
Coefficient of K^2 is (1/4) * SS^4 * (E^2 - 40*SS)
Coefficient of K^1 is E * SS^5
Coefficient of K^0 is SS^6
so after changing the sign of each coefficient we obtain the values from $(2)$.
Code #4: If we know $K$, can we "easily" obtain $x,y,z$?
In other words, are there any "simple relations" among the variables, that would lead to
a quick determination of $x$, when $a,b,c,K$ are known?
Here is the elimination of variables, all but $a,b,c;x,K$.
R.<x,y,z, a,b,c, K,K_inv> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JxK = J.elimination_ideal([y, z, K_inv])
print("Generator(s) of the elimination ideal JxK after eliminating y, z:")
for g in JxK.groebner_basis():
print(f'{g}\n')
And among the many relations shown there is also the following one:
x^3*a^4 - 2*x^3*a^2*b^2 + x^3*b^4 - 2*x^3*a^2*c^2 - 2*x^3*b^2*c^2 + x^3*c^4
- 1/2*x^2*a^4 + x^2*a^2*b^2 - 1/2*x^2*b^4 + x^2*a^2*c^2 + x^2*b^2*c^2 - 1/2*x^2*c^4
+ 4*x^2*a^2*K - 4*x^2*b^2*K - 4*x^2*c^2*K
- 32*x*K^2
+ 8*K^2
This leads to the relation $(3)$.
We can also try to eliminate $K,z$ and obtain relations among $y,z$.
Here are such relations:
sage: for g in J.elimination_ideal([z, K_inv, a, b]).groebner_basis(): print(g, '\n')
x^2*y^2*c^2 + 4*x^2*y*K + 4*x*y^2*K - x^2*K - 2*x*y*K - y^2*K
sage: for g in J.elimination_ideal([z, K_inv, c]).groebner_basis(): print(g, '\n')
x^2*y^2*a^2 + 2*x*y^3*a^2 - 2*x^3*y*b^2 - x^2*y^2*b^2 - 2*x*y^2*a^2 + 2*x^2*y*b^2
+ 4*x^2*y*K - 4*x*y^2*K - x^2*K + y^2*K + 2*x*K - 2*y*K
y^4*a^2 + 2*x^3*y*b^2 + x^2*y^2*b^2 - 2*y^3*a^2 - 2*x^2*y*b^2 + y^2*a^2
+ 8*x*y^2*K - 8*x*y*K - 5*y^2*K + 6*y*K - K
x^4*b^2 + 2*x^3*y*b^2 + x^2*y^2*b^2 - 2*x^3*b^2 - 2*x^2*y*b^2 + x^2*b^2
+ 4*x^2*y*K + 4*x*y^2*K - 4*x^2*K - 8*x*y*K - y^2*K + 4*x*K + 2*y*K - K
Code #5:
Let us see what happens in the case of an isosceles triangle with sides $a=61$, $b=61$, $c=120$.
Its height is $11$, since $61^2-60^2=9^2$.
a, b, c = 61, 61, 120
R.<x,y,z, K,K_inv> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
points = J.variety(ring=AA)
print(f'{a} {b} {c}')
for dic in points:
print(f'x = {dic[x]} y = {dic[y]} z = {dic[z]} K = {dic[K]}')
g = J.elimination_ideal([x, y, z, K_inv]).groebner_basis()[0]
print(f'K is a root of the following polynomial:\n{g.factor()}')
And we obtain:
61 61 120
x = -1.749464267458191? y = -1.749464267458191? z = 4.498928534916382? K = 2449.083313436469?
x = -0.5905364374782210? y = -0.5905364374782210? z = 2.181072874956442? K = 575.6065451903618?
x = 0.1209924404736021? y = 0.1209924404736021? z = 0.7580151190527959? K = 69.52501740622754?
K is a root of the following polynomial:
(1/1800964) * (14884*K^2 + 1757041*K + 52707600) * (121*K^3 - 374400*K^2 + 196020000*K - 11859210000)
A note on the degree of the above polynomial. It is five, not seven.
However, using the formulas for the coefficients in this special case,
a, b, c = 61, 61, 120
s, p, E = (a + b + c)/2, a*b*c, a^2 + b^2 + c^2
SS = s*(s - a)*(s - b)*(s - c)
a7 = 64 * SS * p^2
a6 = -p^4 + 12*p^2*SS*E + 12*SS^2*E^2 + 192*SS^3
a5 = 2 * SS^2 * (E^3 + 13*p^2 + 34*SS*E)
a4 = SS^2 * (p^2*E + 10*SS*E^2 + 103*SS^2)
a3 = SS^3 * (2*p^2 + 17*SS*E)
a2 = (-1/4) * SS^4 * (E^2 - 40*SS)
a1 = - E * SS^5
a0 = - SS^6
var('K')
PI = a7*K^7 + a6*K^6 + a5*K^5 + a4*K^4 + a3*K^3 + a2*K^2 + a1*K + a0
print(PI.factor())
print(PI.roots(ring=AA, multiplicities=False))
we obtain
207360000 * (121*K^3 - 374400*K^2 + 196020000*K - 11859210000) * (14884*K^2 + 1757041*K + 52707600)^2
[69.52501740622754?, 575.6065451903618?, 2449.083313436469?]
So the quadratic factor appears squared.
Code #6:
Let us implement the recursion described in the theoretical section.
def recursion(a, b, c, bit_precision=40):
IR = RealField(bit_precision) # IR is "real field" - precision is given as argument
a, b, c = IR(a), IR(b), IR(c) # pass to the numerical world with the given sides a, b, c
def d2(P, Q):
x1, y1, z1 = P
x2, y2, z2 = Q
x, y, z = x1 - x2, y1 - y2, z1 - z2
return - a^2*y*z - b^2*z*x - c^2*x*y
P = vector(IR, 3, [ IR(1/3), IR(1/3), IR(1/3) ]) # start recursion in centroid
for k in range(10):
print(P)
x, y, z = P
D = vector(IR, 3, [ 0, y/(y+z), z/(y+z)])
E = vector(IR, 3, [ x/(x+z), 0, z/(x+z)])
F = vector(IR, 3, [ x/(x+y), y/(x+y), 0])
PD2, PE2, PF2 = d2(P, D), d2(P, E), d2(P, F)
d2sum = PD2 + PE2 + PF2
X = x * (PE2 + PF2) / 2. / d2sum
Y = y * (PF2 + PD2) / 2. / d2sum
Z = z * (PD2 + PE2) / 2. / d2sum
X, Y, Z = X/(X+Y+Z), Y/(X+Y+Z), Z/(X+Y+Z)
P = vector(IR, 3, [X, Y, Z])
Then calling recursion(6, 9, 13) we return the following result:
sage: recursion(6, 9, 13)
(0.33333333333, 0.33333333333, 0.33333333333)
(0.22960372960, 0.30827505828, 0.46212121212)
(0.20266594188, 0.26478637847, 0.53254767965)
(0.19806794747, 0.26656776311, 0.53536428942)
(0.19727171883, 0.26670350678, 0.53602477439)
(0.19715195703, 0.26675849714, 0.53608954584)
(0.19712950067, 0.26676259568, 0.53610790365)
(0.19712606659, 0.26676417869, 0.53610975472)
(0.19712542143, 0.26676429665, 0.53611028192)
(0.19712532272, 0.26676434216, 0.53611033512)
This "seems to converge" to the algebraic exact value found by code #1.
