Filter löschen
Filter löschen

Is this a bug in eig?

4 Ansichten (letzte 30 Tage)
Kenneth Johnson
Kenneth Johnson am 30 Mär. 2021
Kommentiert: Kenneth Johnson am 31 Mär. 2021
This issue came up with the roots function (see roots_ on File Exchange), where eig was failing on the companion matrix. Following are two test cases; the first one works and the second one doesn't. (I reported this to tech support, but it's not listed in the official Matlab bug reports.)
a = [-1,-1,0;1,0,0;0,1,0];
a(1,3) = 10e-32;
[v,d] = eig(a);
disp(abs(a*v-v*d))
% 1.0e-15 *
% 0.055511151231258 0.055511151231258 0.000000000000000
% 0.166533453693773 0.166533453693773 0.000000000000000
% 0.138777878078145 0.138777878078145 0
disp(num2str(diag(d)))
% -0.5+0.86603i
% -0.5-0.86603i
% 1e-31+0i
a(1,3) = 9e-32;
[v,d] = eig(a);
disp(abs(a*v-v*d))
% 0.000000000000000 0.000000000000000 0.000000000000000
% 0.000000000000000 0.000000000000000 0.000000000000000
% 0.707106781186548 0.707106781186548 0.000000000000000
disp(num2str(diag(d)))
% -0.5+0.86603i
% -0.5-0.86603i
% 0+0i

Akzeptierte Antwort

Matt J
Matt J am 30 Mär. 2021
Bearbeitet: Matt J am 30 Mär. 2021
It's fine. You just need to disable balancing:
a = [-1,-1,0;1,0,0;0,1,0];
for e=[10e-32, 9e-32]
a(1,3) = e;
[v,d] = eig(a,'nobalance');
Discrepancy = abs(a*v-v*d)
end
Discrepancy = 3×3
1.0e+-15 * 0.3608 0.3608 0.0342 0.3554 0.3554 0.0220 0.2355 0.2355 0.2282
Discrepancy = 3×3
1.0e+-15 * 0.3608 0.3608 0.0342 0.3554 0.3554 0.0220 0.2355 0.2355 0.2282
  13 Kommentare
Bruno Luong
Bruno Luong am 31 Mär. 2021
Bearbeitet: Bruno Luong am 31 Mär. 2021
Google I found this page where there is a paper discuss issue with banancing. Someone claims (last post) he has developped a correct balancing method in Lapack 3.5.0.
Kenneth Johnson
Kenneth Johnson am 31 Mär. 2021
For my test case it doesn't work without balancing, but Matlab's balancing doesn't work. There is a discussion of matrix balancing in Numerical Recipes in C++.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by