MATLAB Answers

I am trying to solve for w using determinants

5 views (last 30 days)
Khaled Elmajed
Khaled Elmajed on 23 Jul 2021
Edited: David Goodmanson on 26 Jul 2021
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?

Answers (2)

John D'Errico
John D'Errico on 23 Jul 2021
Edited: John D'Errico on 23 Jul 2021
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.

David Goodmanson
David Goodmanson on 23 Jul 2021
Edited: David Goodmanson on 26 Jul 2021
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.

Community Treasure Hunt

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

Start Hunting!