solve a forced vibration equation for x and w

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

[Sam Chak]

Hi @Fatemeh Salar

Your equation says

that either

or

.

Since the latter is a trivial solution, then I guess you want to solve ?

[Fatemeh Salar]

Hi @Sam Chak

Matter of fact I'm looking for eigenvectors which are given by the solutions xi of

Therefore I had to solved wi in advance by solving below determinant

Now, I'm in a place where I calculated w but don't know how to determine eigenvecotrs which is the solution of xi in first equation.

[Sam Chak]

Hi @Fatemeh Salar

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]

Dear @Sam Chak

Thank you for your sicere effort, but why you decided to define A = -M/K and therefore calculate eig(A) ?

Answer 1

[Sam Chak]

Hi @Fatemeh Salar

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().

[Fatemeh Salar]

@Sam Chak Thank you so much ^^

[Sam Chak]

You are welcome, @Fatemeh Salar. Have a nice day. 😉

Answer 2

[Torsten]

https://de.mathworks.com/help/symbolic/eig.html

First multiply by

 M^(-1) = [1/(3*m) 0; 0 1/m]

Answer 3

[Sam Chak]

Hi @Fatemeh Salar

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)')