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:
matlab precision determinant problem

Subject: matlab precision determinant problem

From: Luka Djigas

Date: 7 May, 2010 14:05:55

Message: 1 of 2

I have the following program

%***************************
format compact; format short g; clear; clc;
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
for i=1:500000
omegan=1.+0.0001*i;
a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0;
a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0;
a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1;
a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2;
if(abs(det(a))<1E-10) sprintf('omegan= %8.3f det= %8.3f',omegan,det(a))
end
end
%***************************

Analytical solution of the above system, and the same program written in fortran
(http://paste-it.net/public/xf94d50/) gives out values of omegan equal to 16.3818 and 32.7636
(fortran values; analytical differ a little, but they're there somewhere).

So, now I'm wondering ... where am I going wrong with this ? Why is matlab not giving the expected
results ?

(this is probably something terribly simple, but it's giving me headaches)

--
Luka

Subject: matlab precision determinant problem

From: Roger Stafford

Date: 7 May, 2010 18:04:28

Message: 2 of 2

Luka Djigas <ldigas@___gmail___.com> wrote in message <3f78u592apubjvsvims31k39akdqka3nrh@4ax.com>...
> I have the following program
>
> %***************************
> format compact; format short g; clear; clc;
> L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
> for i=1:500000
> omegan=1.+0.0001*i;
> a(1,1) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(1,2) = 2; a(1,3) = 0; a(1,4) = 0;
> a(2,1) = 1; a(2,2) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(2,3) = 1; a(2,4) = 0;
> a(3,1) = 0; a(3,2) = 1; a(3,3) = ((omegan^2)*(Jm/(G*J))*d^2)-2; a(3,4) = 1;
> a(4,1) = 0; a(4,2) = 0; a(4,3) = 2; a(4,4) = ((omegan^2)*(Jm/(G*J))*d^2)-2;
> if(abs(det(a))<1E-10) sprintf('omegan= %8.3f det= %8.3f',omegan,det(a))
> end
> end
> %***************************
>
> Analytical solution of the above system, and the same program written in fortran
> (http://paste-it.net/public/xf94d50/) gives out values of omegan equal to 16.3818 and 32.7636
> (fortran values; analytical differ a little, but they're there somewhere).
>
> So, now I'm wondering ... where am I going wrong with this ? Why is matlab not giving the expected
> results ?
>
> (this is probably something terribly simple, but it's giving me headaches)
>
> Luka

  Excuse me for saying so, but that is a very poor test for accuracy, Luka. Why haven't you given 'omegan' the precise values that should create a zero value in the determinant and see how close to zero it actually comes? (Better still, why not put in the theoretical eigenvalues, -2,-1,1,2, in place of that awful expression with G, J, etc. On my machine it gets a precise zero in each case if you do that.)

  According to my computations the following four values (rounded to fifteen places) should theoretically produce a zero determinant:

 omegan = 0.000000000000000
 omegan = 16.381862247021893
 omegan = 28.374217734436371
 omegan = 32.763724494043785

These correspond to indices i of

 i = -10000.000000000000000
 i = 153818.622470218921080
 i = 273742.177344363706652
 i = 317637.244940437842160

(Yes you missed the value omegan = 28.374217734436371 .)

  Whether matlab didn't get the same answer as the fortran program could be as much due to round off errors in fortran as in matlab. The programs were only using inexact integer values for the above i values, so whether they passed the abs(det(a))<1E-10 test was highly dependent on chance roundings in the immediate neighborhood of the 1E-10 threshold value. Matlab might have been right and fortran wrong. Only a test with the symbolic toolbox would determine that.

Roger Stafford

Tags for this Thread

No tags are associated with 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