Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
An interesting numerical problem

Subject: An interesting numerical problem

From: Yi Cao

Date: 22 Oct, 2009 17:27:22

Message: 1 of 5

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

Subject: An interesting numerical problem

From: Matt Fetterman

Date: 22 Oct, 2009 17:48:19

Message: 2 of 5

"Yi Cao" <y.cao@cranfield.ac.uk> wrote in message <hbq4lq$4m1$1@fred.mathworks.com>...
> 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

What happens with N less than 19? If it suddenly doesn't work at 19, that does suggest the memory size? Regards! Matt

Subject: An interesting numerical problem

From: Bruno Luong

Date: 22 Oct, 2009 17:59:07

Message: 3 of 5

"Yi Cao" <y.cao@cranfield.ac.uk> wrote in message <hbq4lq$4m1$1@fred.mathworks.com>...
> 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?
>

cond(P19*P19') is larger than 10^20. Eigen-space are highly parallel. The fact that P19 contains integer does not matter IMO because DET calculation using LU factorization.

>> det(P19*P19')

ans =

    -8

>> det(fliplr(P19*P19'))

ans =

    -1

>> det(flipud(P19*P19'))

ans =

     4

>> det(fliplr(flipud(P19*P19')))

ans =

     1

>>

Subject: An interesting numerical problem

From: Yi Cao

Date: 22 Oct, 2009 19:18:21

Message: 4 of 5

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hbq6hb$5me$1@fred.mathworks.com>...
>
> cond(P19*P19') is larger than 10^20. Eigen-space are highly parallel. The fact that P19 contains integer does not matter IMO because DET calculation using LU factorization.
>

For L*U = P*P', theoretically, L should be P and U should be P' because P is low triangular and P' is upper triangular. Clearly, LU factorization in Matlab is not so smart.

Yi

Subject: An interesting numerical problem

From: someone

Date: 22 Oct, 2009 19:35:23

Message: 5 of 5

"Yi Cao" <y.cao@cranfield.ac.uk> wrote in message <hbqb5s$9st$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hbq6hb$5me$1@fred.mathworks.com>...
> >
> > cond(P19*P19') is larger than 10^20. Eigen-space are highly parallel. The fact that P19 contains integer does not matter IMO because DET calculation using LU factorization.
> >
>
> For L*U = P*P', theoretically, L should be P and U should be P' because P is low triangular and P' is upper triangular. Clearly, LU factorization in Matlab is not so smart.
>
> Yi

IMO, its not so much that "Matlab is not so smart" as it is that due to
floating point arithmetic, numerical problems and theoretical problems
don't always agree. In fact, you can even find "simple" cases where

(a + b) + c ~= a + (b + c)

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us