How can I solve a polynomial system with 17 equation and 17 unknowns?

Hi, I want to solve a 17 unknown system with 17 equations but cannot find the right syntax. Here is what I use so far, but it seems that Matlab is seeking indefinitely, I probably forgot something in there.

a = [2.38705    ;   0.55879 ;   0.07082];
b = [2.58905    ;   0.40879 ;   0.09282];
c = [2.58330    ;   0.70922 ;   0.11628];
d = [2.32884    ;   0.55409 ;   -0.06968];
e = [2.59084    ;   0.34109 ;   -0.13218];
f = [2.59807    ;   0.74465 ;   -0.15085];
g = [2.74650    ;   0.26511 ;   -0.03613];
h = [2.75163    ;   0.75019 ;   -0.03188];

TCP = [2.59000  ;   0.78644 ;   -0.35000];
SP  = [2.59000  ;   0.88644 ;   -0.01211];
WC  = [2.59000  ;   0.78644 ;   -0.01342];

% point c
ch = h - c;
cf = f - c;
ca = a - c;
cb = b - c;

% point f
fh = h - f;
fd = d - f;
fe = e - f;
fWC = WC - f;

% point h
hg = g - h;
hWC = WC - h;
hSP = SP - h;
hTCP = TCP - h;

% axis
WCSP = SP - WC;
WCTCP = TCP - WC;
TCPSP = SP - TCP;

cWC = WC - c;
fSP = SP - f;

TCPz(1) = -0.45;

syms cx cy cz fx fy fz hx hy hz WCx WCy WCz SPx SPy SPz TCPx TCPy

eqn1 = ((cx - hx)^2 + (cy - hy)^2 + (cz - hz)^2 ) == (norm(ch)^2);
eqn2 = ((cx - fx)^2 + (cy - fy)^2 + (cz - fz)^2  ) == (norm(cf)^2);
eqn3 = ((cx - a(1))^2 + (cy - a(2))^2 + (cz - a(3))^2) == (norm(ca)^2);
eqn4 = ((cx - b(1))^2 + (cy - b(2))^2 + (cz - b(3))^2 )== (norm(cb)^2);
eqn5 = ((cx - WCx)^2 + (cy - WCy)^2 + (cz - WCz)^2) == (norm(cWC)^2);

eqn6 = ((fx - hx)^2 + (fy - hy)^2 + (fz - hz)^2 ) == (norm(fh)^2);
eqn7 = ((fx - d(1))^2 + (fy - d(2))^2 + (fz - d(3))^2 ) == (norm(fd)^2);
eqn8 = ((fx - e(1))^2 + (fy - e(2))^2 + (fz - e(3))^2 ) == (norm(fe)^2);
eqn9 = ((fx - WCx)^2 + (fy - WCy)^2 + (fz - WCz)^2 ) == (norm(fWC)^2);
eqn10 = ((fx - SPx)^2 + (fy - SPy)^2 + (fz - SPz)^2 ) == (norm(fSP)^2);

eqn11 = ((hx - g(1))^2 + (hy - g(2))^2 + (hz - g(3))^2 ) == (norm(hg)^2);
eqn12 = ((hx - WCx)^2 + (hy - WCy)^2 + (hz - WCz)^2 ) == (norm(hWC)^2);
eqn13 = ((hx - SPx)^2 + (hy - SPy)^2 + (hz - SPz)^2 ) == (norm(hSP)^2);
eqn14 = ((hx - TCPx)^2 + (hy - TCPy)^2 + (hz - TCPz(1))^2 ) == (norm(hTCP)^2);

eqn15 = ((WCx - SPx)^2 + (WCy - SPy)^2 + (WCz - SPz)^2 ) == (norm(WCSP)^2);
eqn16 = ((WCx - TCPx)^2 + (WCy - TCPy)^2 + (WCz - TCPz(1))^2 ) == (norm(WCTCP)^2);
eqn17 = ((TCPx - SPx)^2 + (TCPy - SPy)^2 + (TCPz(1) - SPz)^2 ) == (norm(TCPSP)^2);

A = solve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12,eqn13,eqn14,eqn15,eqn16,eqn17]);

What am I doping wrong please?

[Walter Roberson]

Could you describe what you are trying to do? It looks like you are trying to find points that have some particular 3D Euclidean distance relationship to each other?

Answer 1

[John D'Errico]

What are you doing wrong?

You are trying to solve something that has been proven mathematically impossible to solve. You have a system of polynomial equations, that are in effect, equivalent to a high order (very high) polynomial, with symbolic coefficients.

It was proved long ago that no solution exists in general if the polynomial is order higher than 4. You problem is easily beyond that point.

Worse, even if you could reduce this to a linear symbolic system, it is way too large to solve, creating millions of terms.

Sorry, but it is easy to create a problem with no solution.

I've not looked through your code, but if you can replace all constants with known values, then you could try to use a numerical solver. But recognize there will probably be hundreds of solutions, and a numerical solver can only find solutions one at a time, depending on your starting values. So you would need to choose zillions of start points, then take the corresponding solutions and probably use clustering to reduce the solution set into a finite number.

[Max_4137]

OK, that would makes sense then. I thought that if you had the same number of unknowns than equations, that would mean that it is solvable.

I will try to review my system then. Thanks for your help.

[John D'Errico]

If it is entirely constant coefficient, then in theory, it has a numerical solution, MAYBE. Or many of them. But finding a symbolic solution can still take a huge amount of effort.

Worse, if the problem reduces to lets say a 5th degree polynomial or higher, even with constant coefficients, there is no analytical solution. Solve is smart enough then to call vpasolve, which will handle the numerical case. But just getting to that point may take a lot of time.

[Max_4137]

Actually when I let it running it will find numerical solutions at the end, just take ages. Is there a way to define some convergence criteria? What I mean is that I have 3 decimals figures, but the solutions it finds are 20 figures after the decimal. I'd like the convergence criteria to be 1e-3 for example or less to make calculations a lot faster. Thanks for your help