why isn't determinant function working?

12 views (last 30 days)
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,:);

Accepted Answer

Cameron Sprent
Cameron Sprent on 16 Feb 2021
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!

More Answers (1)

John D'Errico
John D'Errico on 13 Feb 2021
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
  3 Comments
John D'Errico
John D'Errico on 14 Feb 2021
  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
Cameron Sprent on 14 Feb 2021
  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.

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!