-2

I have reduced my problem to this system of 3 quadratic equations:

  (k1*x + k2*y + k4)^2 = k7      
  (k2*y + k3*z + k5)^2 = k8      
  (k3*z + k1*x + k6)^2 = k9

  (Solve for x, y, z, k1..k9 are constants)

I have spent quite some time not finding an analytical solution. Can someone give me the solution?

P.S. There are some relations between the constants that may simplify the problem. I can post them if needed tomorrow.

Thx Peter

Peter Pohl
  • 17
  • 1
  • 1
    Welcome to StackOverflow. How is this a programming problem? If you add more context, including evidence of your work so far on the problem and a statement of just where you are stuck, this problem would be a better fit for the [Mathematics Stack Exchange](http://math.stackexchange.com/) site. – Rory Daulton Jul 03 '19 at 16:22
  • 1
    I'm voting to close this question as off-topic because it is not about practical computer programming but rather belongs on [Mathematics Stack Exchange](https://math.stackexchange.com/) once the questioner has added more of his own work and explained just where he is stuck. – Rory Daulton Jul 03 '19 at 16:37
  • 1
    I agree that it should be on Math Exchange. However, here's the solution: Take the square root of k7, k8, and k9 to obtain a system of linear equations. Since the square root can be positive or negative, you need to solve the system for all 8 combinations of signs (±sqrt(k7), ±sqrt(k8), ±sqrt(k9)). You will obtain up to 8 real-valued answers or exactly 8 answers if you accept complex coefficients. – Gilles-Philippe Paillé Jul 03 '19 at 18:27

1 Answers1

1

pure math but solvable by programig,... first use the Gilles-Philippe Paillé comment and separate variables:

(k1*x + k2*y + k4)^2 = k7
(k2*y + k3*z + k5)^2 = k8
(k3*z + k1*x + k6)^2 = k9
-------------------------
k1*x + k2*y + k4 = sqrt(k7)
k2*y + k3*z + k5 = sqrt(k8)
k3*z + k1*x + k6 = sqrt(k9)
-------------------------
k1*x + k2*y = sqrt(k7) - k4
k2*y + k3*z = sqrt(k8) - k5
k3*z + k1*x = sqrt(k9) - k6
-------------------------
k1*x + k2*y +  0*z = sqrt(k7) - k4
 0*x + k2*y + k3*z = sqrt(k8) - k5
k1*x +  0*y + k3*z = sqrt(k9) - k6
-------------------------

Now you can rewrite to matrix form

    | k1 k2  0 |
A = |  0 k2 k3 |
    | k1  0 k3 |

    | x |
B = | y |
    | z |

    | sqrt(k7) - k4 |
C = | sqrt(k8) - k5 |
    | sqrt(k9) - k6 |


             A * B =              C
Inverse(A) * A * B = Inverse(A) * C
                 B = Inverse(A) * C

So its simple Inverse of 3x3 matrix. If you expand it to 4x4 by zeor padding and adding 1 to diagonal You can use 4x4 matrix inverse like this:

Just look for matrix_inv in the C++ code example. You can find also the multiplication of matrix and vector there too matrix_mul_vector...

The code would in C++ look like this:

double A[16]=
 {
 k1, 0,k1, 0,
 k2,k2, 0, 0,
  0,k3,k3, 0,
  0, 0, 0, 1
 };
double B[4],C[4]=
 {
 sqrt(k7) - k4,
 sqrt(k8) - k5,
 sqrt(k9) - k6
 };

matrix_inv(A,A);
matrix_mul(B,A,C);

now B should hold your resulting x,y,z values if your equation has solution. What is left is just to add sign combinations as the sqrt of the system loss it ... In case all constants and variables are non negative you can forget about this and use the result directly without trying the 8 combinations ...

If I see it right the combinations are done like this

double C[4]=
 {
 (+/-)sqrt(k7) - k4,
 (+/-)sqrt(k8) - k5,
 (+/-)sqrt(k9) - k6
 };

so for each of the 8 C combinations compute the result ... The combinations itself can be done by a for cycle using 3 lowest bits of iterator variable to decide the sign like:

matrix_inv(A,A);
for (int i=0;i<8;i++)
 {
 if ((i&1)==0) C[0]=+sqrt(k7)-k4; else C[0]=-sqrt(k7)-k4;
 if ((i&2)==0) C[1]=+sqrt(k8)-k5; else C[1]=-sqrt(k8)-k5;
 if ((i&4)==0) C[2]=+sqrt(k9)-k6; else C[2]=-sqrt(k9)-k6;
 matrix_mul(B,A,C);
 // here B holds the i-th solution
 }

In case of complex domain just change the double with a complex data-type...

Spektre
  • 49,595
  • 11
  • 110
  • 380
  • 2
    Thx alot, i did not take the squareroot and got stuck with many ±sqrt-terms when solving one equation for x. I did not see that you could do the ± permutations as last step. Thx to Gilles-Philippe Paillé too. – Peter Pohl Jul 04 '19 at 08:23
  • @PeterPohl its not the last step as you multiply by inverse matrix of A from left after the combination of `+/-` but I got your point :) – Spektre Jul 04 '19 at 14:35
  • 1
    @Spektre Thanks for the detailed answer! I noted a small mistake in the line `Inverse(A) * B = Inverse(A) * C`. It should be `Inverse(A) * A * B = Inverse(A) * C`. Also, I think that using the 4x4 inverse is a bit overkill for this task. I understand that it is to reuse existing code, but the performance can be greatly affected. – Gilles-Philippe Paillé Jul 04 '19 at 15:48
  • @Gilles-PhilippePaillé Nice catch thx I repaired it. Yes 4x4 is slower but not that much. I have also 3x3 somewhere but too lazy to search for it and that 4x4 is already posted. btw I use this [full pseudo inverse matrix](https://stackoverflow.com/a/42293434/2521214) for orthonormal matrices which is much much faster but not usable in this case ... as the ortho-normality is not granted. Here are [GEM and Determinant based inverses](https://stackoverflow.com/a/55439139/2521214) for general case and I also got N-Dimensional matrix/vector template there ... – Spektre Jul 04 '19 at 18:02
  • @Gilles-PhilippePaillé btw that last template was a written as a core math for mine 4D engine related to this QA: [how should i handle (morphing) 4D objects in opengl?](https://stackoverflow.com/a/44970550/2521214) – Spektre Jul 04 '19 at 18:13