Azzera filtri
Azzera filtri

Symbolic solve throwing division by zero error

23 visualizzazioni (ultimi 30 giorni)
The following code periodically throws a division by zero error.
m = 4;
n = 4;
X = sym('X', [m n]);
A = randi([0 1],n,n^2);
B = randi([0 1],m,m^2);
E = X*A - B*kron(X,X);
S = solve(E == 0,X);
S
S = struct with fields:
X1_1: 0 X2_1: 0 X3_1: 0 X4_1: 0 X1_2: 0 X2_2: 0 X3_2: 0 X4_2: 0 X1_3: 0 X2_3: 0 X3_3: 0 X4_3: 0 X1_4: 0 X2_4: 0 X3_4: 0 X4_4: 0
Is it because there are too many variables?
Is it trying to numerically get solutions and is dividing by small numbers?
The error message is an image because the error rarely happens and I could not get it to happen while writing up this question.
  8 Commenti
Joseph Daynger Ruiz
Joseph Daynger Ruiz il 3 Mar 2023
Modificato: Joseph Daynger Ruiz il 3 Mar 2023
The places where zero is the only solution seem to preform really well. In this example with Q, the matrix with the division by zero problem above, the idenity matrix and zero should both be solutions. This has really shaken my confidence in the solve functions abilty to solve. I would like to know what exatly it is doing and if the coeficents being doubles has anything to do with it. Should I convert Q to an integer symbolic matrix first? At the moment I would like to know how a algebraic system comes to the concusion to divide by zero.
The following are examples I found after finding that solve may not be behaving in a typical way.
syms a b x;
solve(x*a - b*x*x == 0,x)
ans = 
This first computation does not provide me with the assumtion that it made that b was nonzero. I thought I recalled in the documentation that solutions also come with the assumtions used. In the past I have used solve to find the exact assuptions where a solve failed and then used that to invesigate the problem spots by hand. The solution should have two solutions for a and b nonzero, one solution when ab=0 but not both zero simultaneously, and infinite solutions when a and b are zero.
The code ran over all cases of zero/one matraceis in the m,n being two case. After seeing the problem fail in the four case I went back to see if solve had not actually answered correctly in some cases.
A = zeros(2,4);
A(1,1) = 1
A = 2×4
1 0 0 0 0 0 0 0
X = sym('X', [2 2]);
E = X*A - A*kron(X,X);
S = solve(E == 0,X);
S
S = struct with fields:
X1_1: [2×1 sym] X2_1: [2×1 sym] X1_2: [2×1 sym] X2_2: [2×1 sym]
Here the solutions should be the diagonal two by two matrices with zero or one in the top left corner and any number in the lower right corner. However solve dose not find all the solutions. I would expect it to find the idenity matrix before it finds the matrix with zeros in all places except a one in the top left corner.
I would offer some huristics for how many solutions between two random matracies there are but now I am not confident that solve is doing what I expect.
Joseph Daynger Ruiz
Joseph Daynger Ruiz il 3 Mar 2023
Modificato: Joseph Daynger Ruiz il 3 Mar 2023
The problem still persists after some preprocessing that would typically make a system of polynomial equations easier to solve.
X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
[0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
[0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
[0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
fgb = reshape(E,1,[]); % Reshapes the polynomials to be a row vector
% since that is the required input for the next
% function.
plys = gbasis(fgb); % Finds the Gröbner basis for the row vector of polynomials
size(plys,2)
ans = 36
vip = symvar(plys);
S = solve(plys == 0,vip); % 36 equations in 16 unknowns
Error using mupadengine/feval_internal
Division by zero.

Error in sym/solve (line 293)
sol = eng.feval_internal('solve', eqns, vars, solveOptions);
S
The division by zero error is so strange since the coefficents are all integers.

Accedi per commentare.

Risposte (1)

Muskan
Muskan il 21 Apr 2023
Hi Joseph,
I tried running your code in R2022b and it does not seem to give any error for that.
Thanks
  1 Commento
Joseph Daynger Ruiz
Joseph Daynger Ruiz il 21 Apr 2023
Modificato: Joseph Daynger Ruiz il 21 Apr 2023
X = sym('X', [4 4]);
Q = [[1,0,0,0,0,-1, 0,0,0,0,-1, 0,0, 0,0,-1];...
[0,1,0,0,1, 0, 0,0,0,0, 0,-1,0, 0,1, 0];...
[0,0,1,0,0, 0, 0,1,1,0, 0, 0,0,-1,0, 0];...
[0,0,0,1,0, 0,-1,0,0,1, 0, 0,1, 0,0, 0]];
E = X*Q - Q*kron(X,X);
fgb = reshape(E,1,[]); % Reshapes the polynomials to be a row vector
% since that is the required input for the next
% function.
plys = gbasis(fgb); % Finds the Gröbner basis for the row vector of polynomials
size(plys,2)
ans = 36
vip = symvar(plys);
S = solve(plys == 0,vip); % 36 equations in 16 unknowns
Error using mupadengine/feval_internal
Division by zero.

Error in sym/solve (line 293)
sol = eng.feval_internal('solve', eqns, vars, solveOptions);
S
Hello @Muskan, I am excited to have another set of eyes looking into the problem. The code above does run, but that is because the problem cases are sparce. The other comments will give a better context to the problem.
Here is a summary. The space of pairs of zero one matracies of size m by m^2 and n by n^2 have problems being solved. Most combinations are solved. When the problem is solved the number of solutions is finite. The error when the solver fails to work confuses me. The solver says there is a divison by zero error.
My own thoughts on the problem. Shouldn't a algebraic solver working on integeral polynomials know when it can or cannot divide by something? I would expect the error to have been a time out error instead. My plan this summer is to investigate systems of polynomials that do the same thing since this division by zero error goes beyond this particular problem.
I realize the code I posted in this comment is showing the problem on a matrix outside the scope of zero one matraices as it has negative entries, but the problem does occur in the original context of the problem. I've been meaning to run the code again to search for the example with only zeros and ones but, I have been a bit busy.

Accedi per commentare.

Categorie

Scopri di più su Loops and Conditional Statements in Help Center e File Exchange

Prodotti


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by