solve a forced vibration equation for x and w

20 ビュー (過去 30 日間)
Fatemeh Salar
Fatemeh Salar 2022 年 6 月 20 日
コメント済み: Sam Chak 2022 年 6 月 22 日
Hi,
I am about to solve this forced-viration equation; considering I already calculate w in the first part but don't know how to calculate x
I tried using this method to calculate w and it worked
m=1;
k=5;
f=10;
M=[3*m 0;0 m];
K=[2*k -k;-k k];
F=[0;f];
syms w
my_eigenvalues=0;
my_eigenvalues=det(K-(w^2)*M);
w=solve(my_eigenvalues,w);
w=double (w)
But when I start to solve the mentioned equation to calculate x , it won't work !!!
syms x
my_eigenvectors=(K-(w^2)*M)*x;
x=solve(my_eigenvectors,x);
x=double (x)
How can I solve this equation ? K-(w^2)*M)*x
  4 件のコメント
Sam Chak
Sam Chak 2022 年 6 月 20 日
Is this approach acceptable?
m = 1;
k = 5;
M = [3*m 0; 0 m];
K = [2*k -k; -k k];
A = -M\K
A = 2×2
-3.3333 1.6667 5.0000 -5.0000
[V, D] = eig(A) % columns of eigenvectors in matrix V and Diagonal matrix D of eigenvalues
V = 2×2
0.6089 -0.3983 0.7933 0.9172
D = 2×2
-1.1620 0 0 -7.1713
Fatemeh Salar
Fatemeh Salar 2022 年 6 月 20 日
Thank you for your sicere effort, but why you decided to define A = -M/K and therefore calculate eig(A) ?

サインインしてコメントする。

採用された回答

Sam Chak
Sam Chak 2022 年 6 月 21 日
I guess you want solve the forced response vibration problem using the Modal Analysis approach. I'll skip the theory (or else it gets super lengthy) and show you directly how to solve in MATLAB. Hope you enjoy this mini tutorial.
m = 1;
k = 5;
f = 10;
M = [3*m 0; 0 m] % mass matrix
M = 2×2
3 0 0 1
K = [2*k -k; -k k] % stiffness matrix
K = 2×2
10 -5 -5 5
F = [0; f] % force vector
F = 2×1
0 10
Mr = sqrtm(M) % mass matrix square root
Mr = 2×2
1.7321 0 0 1.0000
Kt = Mr\K/Mr % calculate the mass-normalized stiffness matrix: inv(Mr)*K*inv(Mr)
Kt = 2×2
3.3333 -2.8868 -2.8868 5.0000
% solve eigenvalue problem to put eigenvectors in matrix V, and eigenvalues in diagonal matrix D
[V, D] = eig(Kt)
V = 2×2
-0.7992 -0.6011 -0.6011 0.7992
D = 2×2
1.1620 0 0 7.1713
W = diag(D); % extract eigenvalues from D
w1 = sqrt(W(1)) % modal frequency 1
w1 = 1.0780
w2 = sqrt(W(2)) % modal frequency 2
w2 = 2.6779
Fr = V'/Mr*F % modal force vector from the eigenvectors
Fr = 2×1
-6.0110 7.9917
S = Mr\V % modal transformation matrix (SUPER IMPORTANT!)... If you get this wrong, then ...
S = 2×2
-0.4614 -0.3470 -0.6011 0.7992
x0 = [0; 0] % initial condition (ic) of the position physical coordinates x(0)
x0 = 2×1
0 0
r0 = S\x0 % convert to ic of modal coordinates using transformation x = S*r, r = S\x
r0 = 2×1
0 0
dx0 = [0; 0] % initial condition (ic) of the velocity physical coordinates dx(0)
dx0 = 2×1
0 0
dr0 = S\dx0 % do it similarly
dr0 = 2×1
0 0
% successfully decouple the physical equations of vibration into two separate modal equations
syms r1(t) r2(t)
eqn1 = diff(r1,t,2) == - (w1^2)*r1 + Fr(1); % r1" + W(1)*r1 = Fr(1), modal equation 1
eqn2 = diff(r2,t,2) == - (w2^2)*r2 + Fr(2); % r2" + W(2)*r2 = Fr(2), modal equation 2
Dr1 = diff(r1,t);
Dr2 = diff(r2,t);
cond = [r1(0)==r0(1), Dr1(0)==dr0(1), r2(0)==r0(2), Dr2(0)==dr0(2)];
eqns = [eqn1, eqn2];
[r1m(t), r2m(t)] = dsolve(eqns, cond) % solve the modal equations to obtain the modal responses
r1m(t) = 
r2m(t) = 
x1 = S(1,1)*r1m + S(1,2)*r2m % Transforming back into the physical coordinates, x = S*r (Superposition principle)
x1(t) = 
x2 = S(2,1)*r1m + S(2,2)*r2m
x2(t) = 
fplot(x1, [0 20], 'r')
hold on
fplot(x2, [0 20], 'b')
hold off
grid on
xlabel('t')
ylabel('\bf{x}')
title('Forced Response')
legend('x_{1}(t)', 'x_{2}(t)')
Conclusion:
The system responses using the Modal Analysis approach is the same as the results from the numerical solution using ode45().
  2 件のコメント
Fatemeh Salar
Fatemeh Salar 2022 年 6 月 21 日
@Sam Chak Thank you so much ^^
Sam Chak
Sam Chak 2022 年 6 月 22 日
You are welcome, @Fatemeh Salar. Have a nice day. 😉

サインインしてコメントする。

その他の回答 (2 件)

Torsten
Torsten 2022 年 6 月 20 日
First multiply by
M^(-1) = [1/(3*m) 0; 0 1/m]

Sam Chak
Sam Chak 2022 年 6 月 20 日
Maybe like this?
function v = dxdt(t, x)
m = 1;
k = 5;
f = 10;
M = [3*m 0; 0 m]; % mass matrix
K = [2*k -k; -k k]; % stiffness matrix
F = [0; f]; % force input matrix
n = size(M, 1);
A = [zeros(n) eye(n); -M\K zeros(n)]; % state matrix
b = M\F;
B = [zeros(n, 1); b]; % input matrix
v = A*x + B; % state-space
end
xo = [0; 0; 0; 0]; % initial condition
tspan = [0 20]; % simulation time
[t, x] = ode45(@dxdt, tspan, xo);
plot(t, [x(:,1) x(:,2)], 'linewidth', 1.5)
grid on
xlabel('t')
ylabel('\bf{x}')
title('Forced Response')
legend('x_{1}(t)', 'x_{2}(t)')

カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by