matlab precision determinant problem

Question

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 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)

Answer 1

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.

Answer 2

New answer:

You can investigate this problem using symbolic equations, which gives me the correct answers:

>> clear all             %# Clear all existing variables
>> format long           %# Display more digits of precision
>> syms Jm d omegan G J  %# Your symbolic variables
>> a = ((Jm*(d*omegan)^2)/(G*J)-2).*eye(4)+...  %# Create the matrix a
       diag([2 1 1],1)+...
       diag([1 1 2],-1);
>> solns = solve(det(a),'omegan')  %# Solve for where the determinant is 0

solns =

                                0
                                0
            (G*J*Jm)^(1/2)/(Jm*d)
           -(G*J*Jm)^(1/2)/(Jm*d)
       -(2*(G*J*Jm)^(1/2))/(Jm*d)
        (2*(G*J*Jm)^(1/2))/(Jm*d)
  (3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d)
 -(3^(1/2)*(G*J*Jm)^(1/2))/(Jm*d)

>> solns = subs(solns,{G,J,Jm,d},{8e7,77,10540,140/3})  %# Substitute values

solns =

                   0
                   0
  16.381862247021893
 -16.381862247021893
 -32.763724494043785
  32.763724494043785
  28.374217734436371
 -28.374217734436371

I think you either just weren't choosing values in your loop close enough to the solutions for omegan or your threshold for how close the determinant is to zero is too strict. When I plug in the given values to a, along with omegan = 16.3819 (which is the closest value to one solution your loop produces), I get this:

>> det(subs(a,{omegan,G,J,Jm,d},{16.3819,8e7,77,10540,140/3}))

ans =

    2.765476845475786e-005

Which is still larger in absolute amplitude than 1e-10.

Answer 3

I put this as an answer because I cannot paste this into a comment: Here's how Matlab calculates the determinant. I assume the rounding errors come from calculating the product of multiple diagonal elements in U.

Algorithm

The determinant is computed from the triangular factors obtained by Gaussian elimination

[L,U] = lu(A) s =  det(L)        
%# This is always +1 or -1  
det(A) = s*prod(diag(U))