49 views (last 30 days)

Show older comments

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,:);

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

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

Start Hunting!