solve a forced vibration equation for x and w
15 views (last 30 days)
Show older comments
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 Comments
Sam Chak
on 20 Jun 2022
Is this approach acceptable?
m = 1;
k = 5;
M = [3*m 0; 0 m];
K = [2*k -k; -k k];
A = -M\K
[V, D] = eig(A) % columns of eigenvectors in matrix V and Diagonal matrix D of eigenvalues
Accepted Answer
Sam Chak
on 21 Jun 2022
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
K = [2*k -k; -k k] % stiffness matrix
F = [0; f] % force vector
Mr = sqrtm(M) % mass matrix square root
Kt = Mr\K/Mr % calculate the mass-normalized stiffness matrix: inv(Mr)*K*inv(Mr)
% solve eigenvalue problem to put eigenvectors in matrix V, and eigenvalues in diagonal matrix D
[V, D] = eig(Kt)
W = diag(D); % extract eigenvalues from D
w1 = sqrt(W(1)) % modal frequency 1
w2 = sqrt(W(2)) % modal frequency 2
Fr = V'/Mr*F % modal force vector from the eigenvectors
S = Mr\V % modal transformation matrix (SUPER IMPORTANT!)... If you get this wrong, then ...
x0 = [0; 0] % initial condition (ic) of the position physical coordinates x(0)
r0 = S\x0 % convert to ic of modal coordinates using transformation x = S*r, r = S\x
dx0 = [0; 0] % initial condition (ic) of the velocity physical coordinates dx(0)
dr0 = S\dx0 % do it similarly
% 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
x1 = S(1,1)*r1m + S(1,2)*r2m % Transforming back into the physical coordinates, x = S*r (Superposition principle)
x2 = S(2,1)*r1m + S(2,2)*r2m
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 Comments
More Answers (2)
Sam Chak
on 20 Jun 2022
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)')
0 Comments
See Also
Categories
Find more on Applications in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!