How to find the integer values of the following system of equations

Given taking integer values from [-10,10], how to find the solution of sysmtem of equations:

1. 
2. 
3. 
4. 
5. 
6. 

[Sam Chak]

I'm just testing out the non-zero integers within the range that satisfy:

Equation (1)

x² + y² + 1 = z²

2² + 2² + 1 = 3²

8² + 4² + 1 = 9²

Answer 1

[John D'Errico]

These are called Diophantine equations, as integer solutions are required. The square terms and products make them nonlinear. As such, they can often be difficult to solve.

Looking at the various equations, we see several things of interest. First, consider the equations:

x1^2 - x2^2 - x3^2 = 1

Swapping terms around, we will see it as:

x1^2 - 1 = x2^2 + x3^2

For ANY value of x1, this reduces to a circle in the (x2,x3) plane, of radius sqrt(x1^2-1). Essentially, it can be viewed as a variation of conic form, actually hyperbolic in nature once we introduce x1.

The same applies to equations 2 and 3. I'm not sure it directly helps us though. At least not immediately.

The domain of the problem is pretty large though, of size

21^9

ans = 7.9428e+11

potential solutions to search. HOWEVER, now we can return to equation 1. As you can see there, we can see the potential solutions can be identified for those three variables. This we can do using brute force, since the problem here is only 3 dimensional. That is...

[X1,X2,X3] = ndgrid(-10:10,-10:10,-10:10);
ind123 = find(X1.^2 - 1 == X2.^2 + X3.^2);

Is this last test valid? I am sure you have learned it is a bad thing to test for exact equality amoug variables in MATLAB. Here however, the variables are ALWAYS integer, and will never exceed 2^53-1. Therefore the test is a valid one. How may solutions did we find?

numel(ind123)

ans = 26

So there are only 26 possible solutions in that domain in the variables {X1,X2,X3}. They are:

X123 = [X1(ind123),X2(ind123),X3(ind123)]

X123 = 26×3

-9 -4 -8 9 -4 -8 -9 4 -8 9 4 -8 -9 -8 -4 9 -8 -4 -9 8 -4 9 8 -4 -3 -2 -2 3 -2 -2

We can solve for all possible solutions to {X4,X5,X6}, using a similar scheme. They will be slightly different of course, because of the -1 instead of the +1 on the right hand side.

[X4,X5,X6] = ndgrid(-10:10,-10:10,-10:10);
ind456 = find(X4.^2 + 1 == X5.^2 + X6.^2);
numel(ind456)

ans = 180

X456 = [X4(ind456),X5(ind456),X6(ind456)];

In this second case, we find 180 possible solutions. Finally, we can recognize that the 2nd and 3rd equations are identical, except for the change in variables. So there are also 180 possible potential solutions in those variables.

X789 = X456;

We have now reduced the search space from 7.9e11 potential solutions, to only

26*180*180

ans = 842400

So it is WAY smaller. In fact, we could now solve this using brute force, thus direct enumeration. I'll do that using a simple approach.

n123 = size(X123,1);
X1 = X123(:,1);
X2 = X123(:,2);
X3 = X123(:,3);

n456 = size(X456,1);
X4 = reshape(X456(:,1),[1,n456]);
X5 = reshape(X456(:,2),[1,n456]);
X6 = reshape(X456(:,3),[1,n456]);

n789 = size(X789,1);
X7 = reshape(X789(:,1),[1,1,n789]);
X8 = reshape(X789(:,2),[1,1,n789]);
X9 = reshape(X789(:,3),[1,1,n789]);

Think carefully about what I did there, as it is important. Now we need to identify the set of solutions that satisfy equations 4,5,6.

indfin = find((X1.*X4 - X2.*X5 - X3.*X6 == 0) & ...
    (X1.*X7 - X2.*X8 - X3.*X9 == 0) & ...
    (X4.*X7 - X5.*X8 - X6.*X9 == 0));
numel(indfin)

ans = 208

So there are 208 possible solutions to all 6 equations. We can write them out now.

[sub123,sub456,sub789] = ind2sub([n123,n456,n789],indfin);
X1234567890 = [X123(sub123,:),X456(sub456,:),X789(sub789,:)]

X1234567890 = 208×9

-9 -4 -8 -4 -1 -4 -8 -4 -7 9 4 8 -4 -1 -4 -8 -4 -7 -9 -4 -8 4 1 4 -8 -4 -7 9 4 8 4 1 4 -8 -4 -7 9 -4 -8 4 -1 -4 8 -4 -7 -9 4 8 4 -1 -4 8 -4 -7 9 -4 -8 -4 1 4 8 -4 -7 -9 4 8 -4 1 4 8 -4 -7 -9 4 -8 -4 1 -4 -8 4 -7 9 -4 8 -4 1 -4 -8 4 -7

Easy peasy. Well, it is, as long as you think about what you are doing.

Answer 2

[Walter Roberson]

You can use ndgrid to construct a grid of all possible x1 x2 x3 values. Use those to express the first equation to see if the condition holds for the triple. Record the triplets it holds for.

Do the same thing for each of the other two equations to come up triplets for x4 x5 x6, and x7 x8 x9

Now you can start substituting in values to the 4th equation to filter out combinations that do not work, giving you a reduced set of combinations. Put the reduced set through #5 to get a further reduced set, and then put what remains through #6 to get a final list that satisfies all of the equations.

[Wallace]

The thing is there are too many solutions satisfying the possible x1 x2 x3 values. What i am trying now is the following: but it seems it takes forever.

for x1=-5:5
    for x2=-5:5
        for x3=-5:5
            for x4=-5:5
                for x5=-5:5
                    for x6=-5:5
                        for x7=-5:5
                            for x8=-5:5
                                for x9=-5:5
                                    a=x1;b=x2;c=x3;d=x4;e=x5;f=x6;g=x7;h=x8;i=x9;
                                    f1=a^2-b^2-c^2;
                                    f2=d^2-e^2-f^2;
                                    f3=g^2-h^2-i^2;
                                    f4=a*d-b*e-c*f;
                                    f5=a*g-b*h-c*i;
                                    f6=g*d-e*h-f*i;
                                    if f1==1 && f2==-1 && f3==-1 && f4==0 && f5==0 && f6==0
                                        disp(['a=',num2str(x1),'b=',num2str(x2),'c=',num2str(x3),])
                                        disp(['d=',num2str(x4),' e=',num2str(x5),' f=',num2str(x6),])
                                        disp(['g=',num2str(x7),' h=',num2str(x8),' i=',num2str(x9)])
                                        
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

[Steven Lord]

If you have values for x1, x2, and x3 in your first set of loops that don't satisfy your first condition (x1.^2-x2.^2-x3.^2 == 1) does it make sense to try going through the rest of the loops and all the calculations they contain? Using a simpler example, if k is too large I don't want to run the loop over m at all. Note that for k = 4 and k = 5, the output contains no mention of m.

for k = 1:5
    disp("k is " + k)
    if k > 3
        continue    
    end
    for m = 1:5
        disp("    m is " + m)
    end
end

k is 1

    m is 1
    m is 2
    m is 3
    m is 4
    m is 5

k is 2

    m is 1
    m is 2
    m is 3
    m is 4
    m is 5

k is 3

    m is 1
    m is 2
    m is 3
    m is 4
    m is 5

k is 4
k is 5

[Walter Roberson]

[x_1, x_2, x_3] = ndgrid(-10:10);

idx1 = find(x_1.^2 - x_2.^2 - x_3.^3 == 1);
X1_sel = x_1(idx1); X2_sel = x_2(idx1); X3_sel = x_3(idx1);

idx2 = find(x_1.^2 - x_2.^2 - x_3.^2 == -1);
X4_sel = x_1(idx2); X5_sel = x_2(idx2); X6_sel = x_3(idx2);

X7_sel = X4_sel; X8_sel = X5_sel; X9_sel = X6_sel;

whos

  Name          Size               Bytes  Class     Attributes

  X1_sel       65x1                  520  double              
  X2_sel       65x1                  520  double              
  X3_sel       65x1                  520  double              
  X4_sel      180x1                 1440  double              
  X5_sel      180x1                 1440  double              
  X6_sel      180x1                 1440  double              
  X7_sel      180x1                 1440  double              
  X8_sel      180x1                 1440  double              
  X9_sel      180x1                 1440  double              
  cmdout        1x33                  66  char                
  idx1         65x1                  520  double              
  idx2        180x1                 1440  double              
  x_1          21x21x21            74088  double              
  x_2          21x21x21            74088  double              
  x_3          21x21x21            74088  double              

So you have a 65 x 180 x 180 search space for equations 4 to 6, which is only about 2.1 million, Is that truly too much to do stepwise elimination?

Hint: create a meshgrid or ndgrid of row indices, 65 x 180 x 180, and use the row indices to index the appropriate _sel variables.

[John D'Errico]

@Walter Roberson - you made one minor error in your analysis. You wrote this:

idx1 = find(x_1.^2 - x_2.^2 - x_3.^3 == 1);

It should have been:

idx1 = find(x_1.^2 - x_2.^2 - x_3.^2 == 1);

There are only 26 solutions to this problem in the domain of interest, not 65. With that fix, your suggestion should produce the same 208 possible solutions as I found overall to the Diophantine system.

Answer 3

[Sam Chak]

Hi @Wallace

Positive integer solutions for Equations 1, 2, 3 are given below. You can populate them to include the negative integers. Negative integers are need to check for Equations 4, 5, 6. Then, you can follow the advice from @Walter Roberson. Hopefully with these integers, you can reduce the computation time.

Positive integers that satisfy Equation 1:

Positive integers that satisfy Equation 2: . (Note: 'any' is valid for 0)

Solutions for Equation 3 are the same as Equation 2.