You're looking for too small of determinant values because Matlab is using a different determinant function (or some other reason like something to do with the floating point accuracy involved in the two different methods). I'll show you that Matlab is essentially giving you the correct values and a better way to approach this problem in general.
First, let's take your code and change it slightly.
format compact; format short g; clear; clc;
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
vals = zeros(1,500000);
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;
vals(i) = abs(det(a));
if(vals(i)<1E-10)
sprintf('omegan= %8.3f det= %8.3f',omegan,det(a))
end
end
plot(1.+0.0001*(1:500000),log(vals))
All that I've done really is logged the values of the determinant for all values of omegan and plotted the log of those determinant values as a function of omegan. Here is the plot:
You notice three major dips in the graph. Two coincide with your results of 16.3818 and 32.7636, but there is also an additional one which you were missing (probably because your condition of the determinant being less than 1e-10 was too low even for your Fortran code to pick it up). Therefore, Matlab is also telling you that those are the values of omegan that you were looking for, but because of the determinant was determined in a different manner in Matlab, the values weren't the same - not surprising when dealing with badly conditioned matrices. Also, it probably has to do with Fortran using single precision floats as someone else said. I'm not going to look into why they aren't because I don't want to waste my time on that. Instead, let's look at what you are trying to do and try a different approach.
You, as I'm sure you are aware, are trying to find the eigenvalues of the matrix
a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0 0 2 -2]];
, set them equal to
-omegan^2*(Jm/(G*J)*d^2)
and solve for omegan. This is how I went about it:
format compact; format short g; clear; clc;
L = 140; J = 77; Jm = 10540; G = 0.8*10^8; d = L/3;
C1 = (Jm/(G*J)*d^2);
a = [[-2 2 0 0]; [1 -2 1 0]; [0 1 -2 1]; [0,0,2,-2]];
myeigs = eig(a);
myeigs(abs(myeigs) < eps) = 0.0;
for i=1:4
sprintf('omegan= %8.3f', sqrt(-myeigs(i)/C1))
end
This gives you all four solutions - not just the two that you had found with your Fortran code (though one of them, zero, was outside of your testing range for omegan ). If you want to go about solving this by checking the determinant in Matlab, as you've been trying to do, then you'll have to play with the value that you're checking the absolute value of the determinant to be less than. I got it to work for a value of 1e-4 (it gave 3 solutions: 16.382, 28.374, and 32.764).
Sorry for such a long solution, but hopefully it helps.
Update:
In my first block of code above, I replaced
vals(i) = abs(det(a));
with
[L,U] = lu(a);
s = det(L);
vals(i) = abs(s*prod(diag(U)));
which is the algorithm that det is supposedly using according to the Matlab docs. Now, I am able to use 1E-10 as the condition and it works. So maybe Matlab isn't calculating the determinant exactly as the docs say? This is kind of disturbing.