I am trying to solve for w using determinants

I am not that really good with Matlab, however, I have a task that requires me to solve using determinants and the code works but eventually I'm getting an output of very strange values

My code (part of it) is as follows:

%Material Parameters Aluminum 6061
E=68.9*10^9 %[elastic modulus n/m2]
rho= 2900 %[density kg/m3]
%Dimensions
L=0.75 %[length m]
CS=6.68*10^-4 %[cross sectional area m2]
I=365180   %[second order moment of inertia mm4]
%Precalculation
S=rho*CS
C=3.51 %[from figure]
syms w
%matrices
K=((E*I)/L^3)*[12 (6*L) -12 (6*L); (6*L) (4*L^2) (-6*L) (2*L^2); -12 (6*L) 12 (-6*L); (6*L) (2*L^2) (-6*L) (4*L^2)]
M=((S*L)/420)*[156 22*L 54 -13*L; 22*L 4*L^2 13*L -3*L^2; 54 13*L 156 -22*L; -13*L -3*L^2 -22*L 4*L^2]
%Determinants
Z=0
Z=det(K-(w^2)*M)
w=solve(Z,w)

Now the idea is that the determinant of (K-(w^2)*M) should equal to zero and I'm solving for w (eigen value problem)

and my output is very strange and doesn't fit here (maybe you can try running the code?)

Where am I going wrong or what am I missing?

Answer 1

[John D'Errico]

Why do you think there is a problem. We can shorten the coefficients of Z a bit, just to look at what is happening.

vpa(Z,3)
ans = 
   2.33e-7*w^8 - 9.42e+13*w^6 + 1.73e+33*w^4 - 7.02e+51*w^2

So Z is an 8th degree polynomial in w. Much of the time, there would be no analytical solution. We can alwaus look at the roots using numerical methods of course.

So, what happens with solve? We expect 8 roots.

wsol = solve(Z,w,'maxdegree',8);
   % a lengthy mess of characters....

But, a lengthy mess that vpa can resolve.

vpa(wsol,10)
ans =
                                      0
                                      0
- 19617221100.0 + 0.000000001019512392i
     2441609204.0 - 0.0000003591893407i
     3623418836.0 + 0.0000002475562557i
  19617221100.0 - 0.000000001019512392i
   - 2441609204.0 + 0.0000003591893407i
   - 3623418836.0 - 0.0000002475562557i

There are two real roots at w == 0. All other roots are complex. There are NO non-zero real roots. Of course, the imaginary parts of those roots are pretty tiny in comparison to the real parts. My guess is the way this was constructed, those are all essentially real roots. Looking at the coefficients of the polynomial, that makes complete sense. And then when we look at the polynomial itself, you should realize that we could have written the problem in a related form as a 4th degree polynomial. (I could also have divided by w^2 also to deflate the double root at 0.)

And if we display that in a readable form, though this time with 16 digits in the coefficients...

Zhat = simplify(subs(Z, w, sqrt(u)));
vpa(Zhat,16)
ans =
0.0000002331181919822607*u^4 - 94162494878960.04*u^3 + 1.730905943624151e+33*u^2 - 7.021671728455703e+51*u

which has its roots at:

vpasolve(Zhat)
ans =
                                0
5961455502964345883.6107318719232
13129164062555735516.306228816905
384835363551770922074.98206986208

So the original polynomial had roots at +/-:

>> sqrt(vpasolve(Zhat))
ans =
                                0
2441609203.5713548816199079241373
3623418836.2036944508869271368416
19617221096.571525415065736091503

WTP? Looks ok to me. The imaginary parts from the direct solve are spurious, floating point trash.

Effectively, you did nothing wrong. You just do not know how to look at the solutions, and reduce them to something useful.

Answer 2

[David Goodmanson]

Hi Khaled,

complex results do not seem right in this situation.

M-M'
ans =
     0     0     0     0        % symmetric
     0     0     0     0
     0     0     0     0
     0     0     0     0

K-K'
ans =
            0            0            0            0
            0            0  -5.3677e+17            0
            0   5.3677e+17            0            0
            0            0            0            0

K =
   7.1569e+17   2.6838e+17  -7.1569e+17   2.6838e+17
   2.6838e+17   1.3419e+17  -2.6838e+17   6.7096e+16
  -7.1569e+17   2.6838e+17   7.1569e+17  -2.6838e+17
   2.6838e+17   6.7096e+16  -2.6838e+17   1.3419e+17

K(3,2) and K(2,3) have opposite signs, so K is not symmetric. I have my doubts that that is what was intended. To make K symmetic, either both of those elements are positive or both are negative. Making them both positive did not improve matters, but if both are negative, then

%Material Parameters Aluminum 6061
E=68.9*10^9 %[elastic modulus n/m2]
rho= 2900 %[density kg/m3]
%Dimensions
L=0.75 %[length m]
CS=6.68*10^-4 %[cross sectional area m2]
I=365180   %[second order moment of inertia mm4]
%Precalculation
S=rho*CS
C=3.51 %[from figure]
syms w
%matrices
K=((E*I)/L^3)*[12 (6*L) -12 (6*L); (6*L) (4*L^2) (-6*L) (2*L^2); -12 (6*L) 12 (-6*L); (6*L) (2*L^2) (-6*L) (4*L^2)]

% ----------------------
K(3,2) = -K(3,2);
% ---------------------
M=((S*L)/420)*[156 22*L 54 -13*L; 22*L 4*L^2 13*L -3*L^2; 54 13*L 156 -22*L; -13*L -3*L^2 -22*L 4*L^2]
%Determinants
Z=0
Z=det(K-(w^2)*M);
ww = solve(Z,w);
format short g
ww = double(ww)

ww =
            0
            0
      -26.548
  -1.8569e+10
       26.548
   1.8569e+10
  -5.4365e+09
   5.4365e+09

Since the equation is fourth degree in w^2, there should be pairs of roots that are +- of each other, which there are.

There is at least one more concern, though. Do resonant frequencies in the gigaHz make sense for this system? Looking back, you state that the moment of inertia is in mm^4 and I did not see a later adjustment for that. Converting the parameter to units of m^4 results in

ww =
            0
            0
  -1.7903e-05
       -18569
   1.7903e-05
        18569
      -5436.5
       5436.5

which seems to be more in the ballpark.