# Solve a system of 4-th order polynomial equations (numerically)

I have a function defined on $$S^6 \times S^6$$ (two spheres embedded in $$\mathbb{R}^7$$ individually). Let us call this vector $$\vec{f}(\vec{x},\vec{y})$$, where $$\vec{x}$$ and $$\vec{y}$$ are points in the two spheres. The target space of the function if $$\mathbb{R}^{13}$$, namely $$\vec{f}: S^6 \times S^6 \rightarrow \mathbb{R}^{13}$$. The vectors $$\vec{x}, \vec{y}$$ live in the spheres, so $$\vec{x} \cdot \vec{x} = 1$$ and $$\vec{y} \cdot \vec{y} = 1$$.

This is the problem I am trying to solve. I want to find at least two points $${\vec{x}_0,\vec{y}_0}$$ such that

$$\vec{f}(\vec{x}_0, \vec{y}_0) = \vec{0}$$

Each entry $$\vec{f}$$ is a polynomial of the 4-th order in the entries of the vectors and can be written formally as

$$f^i(\vec{x}, \vec{y}) = \sum_{jknm} a^{i}_{jknm} x^j x^k y^n y^m$$ where $$a^{i}_{jknm}$$s are real coefficients.

I have saved the vector $$\vec{f}$$ in this file .wbx (download it) that you can import with the command

Uncompress @ Import["vector.wbx", "String"]


In this file, the entries of the vectors $$\vec{x}$$ and $$\vec{y}$$ are labelled with $$x_i$$ and $$y_i$$ for $$i = 1, 2, ..., 7$$.

Given the explicit values of the coefficients $$a^i_{jknm}$$, I have tried to solve the system with Mathematica. Obviously, the system is very difficult and an analytic solution is hard to find. So, I have tried to find the solution numerically, but Mathematica takes too long and does not give an answer.

I have tried with

Solve[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x, y} // Flatten]
NSolve[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x , y} // Flatten]


and also I have looked if Mathematica can find just a solution with

FindInstance[{f == Table[0, {i, 1, 13}, x.x == 1, y.y == 1}, {x, y} // Flatten]


but still it takes too long.

Do you have any suggestion in order to solve this problem with Mathematica?

• What about FindRoot? – Henrik Schumacher Nov 10 '18 at 17:32
• @HenrikSchumacher thanks for the suggestion. The point is that I have 13 equations from the system + the 2 constraints $\vec{x}.\vec{x}=1$, $\vec{y}.\vec{y}=1$ which gives 15 equations. On the contrary, I have only 14 variables. So, FindRoot complains that the numbers do not match. – apt45 Nov 10 '18 at 17:43
• @HenrikSchumacher Even If I add a fictitious variable k in FindRoot with arbitrary initial value, it finds a singular Jacobian at the initial condition points. See here i63.tinypic.com/2is6omd.png . I do not why. – apt45 Nov 10 '18 at 17:54
• Possible duplicate: mathematica.stackexchange.com/a/91247/4999 – Michael E2 Nov 10 '18 at 22:17

Perhaps the easiest solution, inferred from a comment by the OP below:

obj = Uncompress@Import["/tmp/vector.wdx", "String"];

obj /. {x3 -> 0, x4 -> 0, x5 -> 0, x6 -> 0, x7 -> 0,
y3 -> 0, y4 -> 0, y5 -> 0, y6 -> 0, y7 -> 0}
(*  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}  *)


So x1, x2, y1, y2 may take on any values.

Original numerical approach

One could minimize the squares of the residuals. Use higher WorkingPrecision for a more accurate result.

obj = Uncompress@Import["/tmp/vector.wdx", "String"];
vars = Variables@obj;
sumOfSq = Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} // #.# &;

SeedRandom (* for reproducibility *)

{acc, sol1} = FindMinimum[
(*
{1.87883*10^-32,
{x1 -> 0.770225, x2 -> 0.637772, x3 -> -1.07917*10^-10, x4 -> -1.40295*10^-9,
x5 -> 1.1682*10^-9, x6 -> 2.92351*10^-8, x7 -> -1.02405*10^-8,
y1 -> 0.626869, y2 -> -0.779124, y3 -> 1.54472*10^-8, y4 -> -5.33267*10^-12,
y5 -> -1.04004*10^-9, y6 -> -5.10446*10^-10, y7 -> -2.39731*10^-10}}
*)

{acc, sol2} = FindMinimum[
WorkingPrecision -> 32, MaxIterations -> 500]
(*
{7.6898897201977925978778380482460*10^-64,
{x1 -> -0.90297104074904814666962995570081,
x2 -> 0.42970140745473574888183421474225,
...,
y1 -> 0.65977018314385699231516512470880,
y2 -> 0.75146743471312275002694856178381,
...,
y7 -> -2.6423411962015994155963564116132*10^-18}}
*)

Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} /. sol1 // Abs // Max
Flatten@{obj, #.# - 1 & /@ Partition[vars, 7]} /. sol2 // Abs // Max
(*
1.09623*10^-16
0.*10^-32
*)


Note that sol1 is pretty accurate at machine precision, just as accurate as sol2 when it is converted to machine precision:

Flatten@{obj, (#.# - 1) & /@ Partition[vars, 7]} /. N@sol2 // Abs // Max
(*  1.38778*10^-16  *)

• Thanks. Can you find another solution? For some reasons I didn't explain, I was expecting to find a solution for all values of x1 x2 y1 y2 and the other variables setted to zero. Can you find a solution for, say, X5 y5 != 0? – apt45 Nov 11 '18 at 2:54
• Do you think this method can be generalized by forcing these additional constraints? – apt45 Nov 11 '18 at 2:54
• @apt45 Every time I run the code with a new random starting vector, I get a different solution with the higher index variables (e.g. x5 etc.) close to zero. Note sol1 and sol2 are different, and for sol2 the higher index variables are much smaller than sol1. It seems the numerics limits the accuracy of these zero value to one-half of working precision, i.e. $10^{-wp/2}$. – Michael E2 Nov 11 '18 at 3:46
• Thank you very much. So, are you saying that you expect the only solution to be for higher index variables very close to zero? – apt45 Nov 11 '18 at 10:51
• @apt45 Yes, if they are computed with floating-point. Even a one-bit rounding error would give a relative error of $dx/x\approx2 \times10^{-16}$. In practice, one should expect more than that, because as errors propagate through a computation, relative error is more likely to increase than decrease. The error in the function value should be around $\|df\|.dx$ at most where $\|df\|$ is the (matrix) norm of the derivative, so one should not expect the zeros to be produce zero exactly. There is a somewhat technical, related discussion here. – Michael E2 Nov 11 '18 at 14:08