Prev: Error in Matlab communication with Com ports
Next: Converting date to serial date (Beginner problem)
From: Luka Djigas on 7 May 2010 10:05
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
From: Roger Stafford on 7 May 2010 14:04
Luka Djigas <ldigas@___gmail___.com> wrote in message <3f78u592apubjvsvims31k39akdqka3nrh(a)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
 | 
Pages: 1
Prev: Error in Matlab communication with Com ports
Next: Converting date to serial date (Beginner problem)