why isn't determinant function working?

I have a matrix for which I know the determinant is equal to 0, and the matrix contains an unkown - omega. I'm trying to find the omega values that make the determinant equal to zero. I know what these values are as they are given in the paper I got the matrix from. Furthermore, mathematica found the correct values (results = [0.599924;1.55627;2.8629;4.42236;6.21874;8.65986;11.0818;13.1635;15.2926]). I have checked the entries over and over and I'm sure there aren't any typos. Why isn't this working? I've also tried to plot the value of the determinant for different values of omega and the plot is not correct. I'm happy to supply that script too if that's helpful.

clear variables

% plate properties
h = 0.01; %plate thickness
a = 0.1; %outer radius
d = 0.02; %inner radius
mu = 0.3; %poisson's ratio

%given properties in paper
beta = d/a;
K_sq = 0.66;  %given in paper - not sure what it is
k_sq = (1-mu)*K_sq/2; %dimensionless shear coefficient

syms omega r MAT1
%set up the functions defined in paper
alpha_sq = (1/12)*(h/a)^2;
delta_1_sq = (omega^2/2)*(1+(1/k_sq)+sqrt(((1-(1/k_sq))^2 + (4/(alpha_sq*omega^2)))));
delta_2_sq = (omega^2/2)*(1+(1/k_sq)-sqrt(((1-(1/k_sq))^2 + (4/(alpha_sq*omega^2)))));
delta_1 = sqrt(delta_1_sq);
delta_2 = sqrt(delta_2_sq);
sigma_1 = delta_2_sq*(omega^2 - k_sq/alpha_sq)^-1;
sigma_2 = delta_1_sq*(omega^2 - k_sq/alpha_sq)^-1;

%Matrix
MAT1(1,1) = besselj(0,delta_1);
MAT1(1,2) = besselj(0,delta_2);
MAT1(1,3) = bessely(0,delta_1);
MAT1(1,4) = bessely(0,delta_2);
MAT1(2,1) = delta_1*(1-sigma_1)*besselj(1,delta_1);
MAT1(2,2) = delta_2*(1-sigma_2)*besselj(1,delta_2);
MAT1(2,3) = delta_1*(1-sigma_1)*bessely(1,delta_1);
MAT1(2,4) = delta_2*(1-sigma_2)*bessely(1,delta_2);
MAT1(3,1) = delta_1*(1-sigma_2)*besselj(1,(delta_1*beta));
MAT1(3,2) = delta_2*(1-sigma_2)*besselj(1,(delta_2*beta));
MAT1(3,3) = delta_1*(1-sigma_1)*bessely(1,(delta_1*beta));
MAT1(3,4) = delta_2*(1-sigma_2)*bessely(1,(delta_2*beta));
MAT1(4,1) = delta_1*besselj(1,(delta_1*beta));
MAT1(4,2) = delta_2*besselj(1,(delta_2*beta));
MAT1(4,3) = delta_1*bessely(1,(delta_1*beta));
MAT1(4,4) = delta_2*bessely(1,(delta_2*beta));

%find the frequencies
i =9; %no. of modes to find
results=zeros(i,1);
DET1 = det(MAT1);

for k = 1:5
    results(k,1) = vpasolve(DET1,omega,[0 20],'Random',true);
end

[~,idx] = sort(results(:,1));
sortedresults = results(idx,:);

Answer 1

[Cameron Sprent]

So I think the errors are due to rounding errors in the det function. When I find the determinant symbolically and only substitute the numbers in right before I try to find the roots, I get the correct roots. Simple solution, just frustrating!

Answer 2

[John D'Errico]

I think your problem may be due to the imaginary part of this rather nasty function.

F = matlabFunction(DET1);

F(1)
ans =
      -27.614 - 9.0262e-13i
      
F(5)
ans =
        40512 + 4.9958e-07i
        
fplot(@(om) imag(F(om)),[0 20])

As you can see, this mess has an imaginary part that is HIGHLY oscillatory, but almost entirely effectively zero. If instead we look at the real part, we see:

fplot(@(om) real(F(om)),[0 20])

So there are now multiple roots in that region, IF we ignore that imaginary part.

vpasolve(real(DET1),omega,[0 20])
ans =
2.8380813015785161825748287968709

[Cameron Sprent]

Thanks for replying so quickly! I have managed to plot the determinant in 3d - input as a real number, and splitting the output into the real and imaginary parts gives a nice way of visualing it. However, the zeros of that plot do not occur at the points where I expect them to. I guess the two questions ive got are a) why is the plot wrong and b) Why is it that mathematica can find the roots but matlab cannot? thank you!

[John D'Errico]

1. How do I know if you implemented the equations properly? Start with that.
2. Verify that if you substitute the value of omega that you THINK is a solution, does it return zero? If not, then there is a probem with the equation you posed. What you may have done wrong, I cannot tell.
3. Next, remember that a solver is trying to find a point where the value is ZERO. And that means both a real and imaginary part equal to zero.

[Cameron Sprent]

1. The equations are correct.
2. When I substitue the correct values of omega it does not return zero.
3. I understand that, but as you can see on the attached plot the points where the determinant is zero are not correct.

Is there potentially an issue with scalability of bessel functions? I saw that you can scale the function, but it doesn't seem to work if the function is used symbolically.