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))<1E10) 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://pasteit.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))<1E10 test was highly dependent on chance roundings in the immediate neighborhood of the 1E10 threshold value. Matlab might have been right and fortran wrong. Only a test with the symbolic toolbox would determine that.
Roger Stafford
