Path: news.mathworks.com!not-for-mail
From: "Yi Cao" <y.cao@cranfield.ac.uk>
Newsgroups: comp.soft-sys.matlab
Subject: An interesting numerical problem
Date: Thu, 22 Oct 2009 17:27:22 +0000 (UTC)
Organization: Cranfield University
Lines: 25
Message-ID: <hbq4lq$4m1$1@fred.mathworks.com>
Reply-To: "Yi Cao" <y.cao@cranfield.ac.uk>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1256232442 4801 172.30.248.38 (22 Oct 2009 17:27:22 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 22 Oct 2009 17:27:22 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 69713
Xref: news.mathworks.com comp.soft-sys.matlab:579391


I have met an interesting numerical problem:

Let's define a function first:

function P=polyconvert(n)
P = eye(n);
for k=1:n
     P(k,1:k) = fliplr(poly(ones(1,k-1)));
end

Then use the function as follows:

P19 = polyconvert(19);

This is a low triangular matrix with all diagonal elements equal to 1. Therefore, the determinant of the matrix should be 1, as well any product of P19 and P19'. However, in Matlab 2009b, I got the following results:

det(P19)        % I got 1, it is fine.
det(P19'*P19) % I also got 1.
det(P19*P19') % I got -8!

It is very strange, how can det(P19*P19') equal to -8? All elements of P19 are integers, hence should not have any error here. The maximum elements of P19*P19' is only about 1e10. Although it is slightly over the 2^32, but that should not be a problem for double floating point calculation, isn't it? Anyone can give a clue what is going wrong here?

Thanks for help.

Yi