Apostila_MCAE (2)

Prévia do material em texto

- APOSTILA - 
(MÉTODOS COMPUTACIONAIS APLICADOS À ENGENHARIA) 
 
 
 
 
 
 
 
 
 
 
 
ATIVIDADE 1 
 
1- Dadas as matrizes A e B: 
 
Onde: 
𝐴 = [
1 2 −2
1 0 3
−4 5 10
], 𝐵 = [
−1 4 −5
10 −3 2
5 2 −1
]; 
 pede-se determinar: 
C = A + B; 
D = A * B; 
E = inversa de A; 
F = Determinante de B; 
 
%Programa Matrizes 
 
a = [1 2 -2;-1 0 2; 4 -5 3]; 
b = [-1 1 -3;2 2 2;-2 4 -1]; 
 
c = a+b 
 
d = a*b 
 
e = inv(a) 
 
f = det(b) 
 
2- Dados os vetores 1 2
2 4
3 2
2 2
 G e G
−   
   
= = −
   
   −   
 
Determinar: 
Os produtos vetoriais: 
 
1 1 2 2 2 1 h G G e h G G=  = 
 
O módulo de 1h 
O produto escalar: 
 
 
1 2 . gg G G= 
 
%Produto Vetorial 
 
g1 = [2;3;-2]; 
 
g2 = [-4;-2;2]; 
 
h1 = cross(g1,g2) 
 
h2 = cross(g2,g1) 
 
modh1 = sqrt(h1(1)^2+h1(2)^2+h1(3)^2) 
 
%Produto Escalar 
gg = g1’*g2 
 
 
3 - Dado o polinômio: 
 
5 4 3( ) 3 - 7 -3 1= p x x x x x+ + 
 
Determinar suas raízes. 
 
Resolver as seguintes equaçoes: 
 
2 22 0 2 0 - x e x+ = − = 
 
 
%Programa - Polinômio 
 
clear all 
 
p = [1 3 -7 0 -3 1]; 
q = roots(p) 
 
p1 = [-1 0 2]; 
q1 = roots(p1) 
 
p2 = [1 0 -2]; 
q2 = roots(p2) 
 
4- Determinar a solução do sistema (1). 
 
{
𝑥 + 𝑦 − 3𝑧 + 𝑤 = 7
−2𝑥 − 𝑦 + 𝑧 + 4𝑤 = 9
𝑥 + 3𝑦 − 𝑧 − 4𝑤 = −4
4𝑥 − 𝑦 + 3𝑧 − 𝑤 = −5
 (1) 
O sistema (1) escrito na forma matricial é dado por (2). 
1 1 3 1 7
2 1 1 4 9
1 3 1 4 4
4 1 3 1 5
x
y
z
w
−     
     
− −
     =
     − − −
     
− − −     
 (2) 
A Eq. (2) é da forma (3). 
[A][X]=[B] (3) 
Pré-multiplicando (3) por [A]-1 , tem-se (4). 
 
[X]= [A]-1 [B] (4) 
 
%Solução do Sistema 
clear all 
 
a = [1 1 -3 1;-2 -1 1 4;1 3 -1 -4;4 -1 3 -1]; 
 
b = [7;9;-4;-5]; 
 
x = inv(a)*b 
ATIVIDADE 2 
 
1 - Dadas as funções: 𝑦1 = 𝑥2 + 4 e 𝑦2 = 𝑥 + 5, elabore um programa 
computacional, para gerar os gráficos de: 𝑦1 em função de 𝑥; 𝑦2 em 
função de 𝑥, com 𝑥𝑖𝑛𝑖𝑐𝑖𝑎𝑙 = −4, com incremento 𝑑𝑥 = 0,1 e um total de 
100 pontos. 
%Gráfico XY 
 
clear all 
close all 
 
x = -4; 
dx = 0.1; 
 
for i = 1:1:100; 
 
 y1 = x^2+4; 
 y2 = x+5; 
 
 
 xv(i) = x; 
 y1v(i)=y1; 
 y2v(i)=y2; 
 x = x+dx; 
end 
 
figure(1) 
plot(xv,y1v,xv,y2v);grid; 
legend('y1','ý2') 
xlabel('x') 
ylabel('y' 
 
 
2 - Dada a função 2z xy= , elabore um programa computacional, para 
gerar o gráfico de z em função de 𝑥 𝑒 𝑑𝑒 𝑦, com xinicial = yinicial = 0 e 
incremento 𝑑𝑥 = 𝑑𝑦 = 0,1 e um total de 200 pontos para x e 100 pontos 
para y. 
%Gráfico XYZ 
 
clear all 
close all 
 
x = 0; 
dx = 0.1; 
y = 0; 
dy = 0.1; 
 
for j = 1:1: 100 
 x = 0; 
 for i = 1:1:200 
 
 z = x*y^2; 
 xv(i,j) = x; 
 yv(i,j) = y; 
 zv(i,j) = z; 
 
 x = x+dx; 
 end 
 y = y + dy; 
end 
-4 -2 0 2 4 6
0
5
10
15
20
25
30
35
40
x
y
 
 
y1
ý2
plot3(xv,yv,zv,'k'); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
 
 3 - Elaborar um programa para geração da função pulso de amplitude de 
10 N, 
 com largura de pulso de 5 segundo. 
 
%Gráfico ForçaxTempo 
 
clear all 
close all 
 
t = 0.0; 
dt = 0.1; 
 
for i = 1:1:500; 
 if i>=0&i<=50 
 u = 10.0; 
 else 
 u=0; %mudar para função degrau (u = 10) 
 end 
 tv(i) = t; 
 uv(i) = u; 
 t = t+dt; 
end 
 
0
5
10
15
20
0
5
10
0
500
1000
1500
2000
 X (m) Y (m)
 Z
 (
m
)
figure(1) 
plot(tv,uv);grid; 
xlabel('Tempo (s)') 
ylabel('Força (N)') 
 
 
 
4 - Elaborar um programa para geração o gráfico dado pela seguinte 
função: 
Se t for maior e igual a zero e menor e igual a 2 seg, u = -5t; 
Se t for maior que 2 e menor e igual a 5 seg, u = -10, e 
Se t for maior que 5 e menor e igual a 8 seg, u = 3. 
%Gráfico ReferênciaxTempo 
 
clear all 
close all 
t = 0.0; 
dt = 0.01; 
 
for i = 1:1:801 
if t>=0 & t<=2 
 u = -5*t; 
 elseif t>2 & t<=5 
 u = -10; 
0 10 20 30 40 50
0
2
4
6
8
10
Tempo (s)
F
o
rç
a
 (
N
)
 
 elseif t>5 & t<=8 
 u = 3; 
 end 
 uv(i)=u; 
 tv(i) = t; 
 t = t+dt; 
end 
plot(tv,uv);grid 
xlabel('Tempo (s)') 
ylabel('Referencia (N)') 
legend('Referencia') 
 
 
 
 
 
 
 
 
 
 
 
0 1 2 3 4 5 6 7 8
-10
-8
-6
-4
-2
0
2
4
Tempo (s)
R
e
fe
re
n
c
ia
 (
N
)
 
 
Referencia
ATIVIDADE 3 
 
1 – Determinar o espaço de trabalho do robô de um grau de liberdade 
(1gdl) quando o motor girar de 0 a 360 graus. 
 
% Espaço de trabalho - 1 gdl (360 graus) 
 
clear all 
close all 
 
rd = pi/180.; 
g = 1/rd; 
 
r = 0.3; 
 
teta = 0*rd; 
dteta = 1*rd; 
 
for i = 1:1:361 
 
 rv = [r*cos(teta);r*sin(teta);0]; 
 x(i) = rv(1); 
 y(i) = rv(2); 
 teta = teta+dteta; 
end 
figure(1) 
plot(x,y,'k') 
xlabel(' X (m)') 
ylabel(' Y (m)') 
legend('espaço de trabalho') 
 
 
2 – Determinar a trajetória do órgão efetuador (partícula A) do robô de 
1gdl, quando o motor girar de 0 a 90 graus. 
% Espaço de trabalho - 1 gdl (90 graus) 
 
clear all 
close all 
 
rd = pi/180.; 
g = 1/rd; 
 
r = 0.3; 
 
teta = 0*rd; 
dteta = 1*rd; 
 
for i = 1:1:91 
 
 rv = [r*cos(teta);r*sin(teta);0]; 
 x(i) = rv(1); 
 y(i) = rv(2); 
 teta = teta+dteta; 
end 
figure(1) 
plot(x,y,'k') 
xlabel(' X (m)') 
-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
 X (m)
 Y
 (
m
)
 
 
espaço de trabalho
ylabel(' Y (m)') 
legend('trajetoria') 
 
 
 
3 – Determinar o espaço de trabalho do robô de dois graus de liberdade 
(2gdl) quando o motor girar de 0 a 360 graus, e deslocamento do raio de 
0,1 a 0,3. 
 
 
-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
 X (m)
 Y
 (
m
)
 
 
trajetoria
% Espaço de trabalho - 2 gdl (360 graus e deslocamento do raio) 
 
clear all 
close all 
 
rd = pi/180.; 
g = 1/rd; 
 
r = 0.1; 
dr = 0.005; 
 
teta = 0*rd; 
dteta = 1*rd; 
 
for j = 1:1:40 
for i = 1:1:361 
 
 rv = [r*cos(teta);r*sin(teta);0]; 
 x(i,j) = rv(1); 
 y(i,j) = rv(2); 
 teta = teta+dteta; 
end 
r = r + dr; 
end 
figure(1) 
plot(x,y,'k') 
xlabel(' X (m)') 
ylabel(' Y (m)') 
legend('espaço de trabalho') 
 
 
4 – Determinar a trajetória do órgão efetuador (partícula A) do robô de 
2gdl, quando o motor girar de 0 à 90 graus e o pistão pneumático se 
deslocar de 0,1 à 0,2 m. 
% Espaço de trabalho - 2 gdl (90 graus e deslocamento do raio) 
 
clear all 
close all 
 
rd = pi/180.; 
g = 1/rd; 
 
r = 0.1; 
dr = 0.00111; 
 
teta = 0*rd; 
dteta = 1*rd; 
 
for i = 1:1:91 
 
 rv = [r*cos(teta);r*sin(teta);0]; 
 x(i) = rv(1); 
 y(i) = rv(2); 
 teta = teta+dteta; 
-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
 X (m)
 Y
 (
m
)
 
 
espaço de trabalho
r = r + dr; 
end 
figure(1) 
plot(x,y,'k');grid 
xlabel(' X (m)') 
ylabel(' Y (m)') 
legend('trajetoria') 
 
5 – Determinar o espaço de trabalho do robô de dois graus de liberdade 
(2gdl), quando o motor girar de 0 a 360 graus e a altura deslocar de 0 à 
0,6 m. 
 
-0.02 0 0.02 0.04 0.06 0.08 0.1 0.12
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
 X (m)
 Y
 (
m
)
 
 
trajetoria
% Espaço de trabalho - 2 gdl (360 graus e deslocamento de altura) 
 
clear all 
close all 
 
r = 0.5; 
 
rd = pi/180.0; 
g = 1/rd; 
 
dteta = 1*rd; 
h = 0.0; 
dh = 0.01; 
 
for j =1:1:60 
 
 teta = 0.0*rd; 
 
 for i = 1:1:360 
 z = h; 
 rav = [r*cos(teta);r*sin(teta);z]; 
 
 ravx(i,j) = rav(1);ravy(i,j) = rav(2); 
 ravz(i,j) = rav(3); 
 
 teta = teta + dteta; 
 end 
 h =h + dh; 
end 
 plot3(ravx,ravy,ravz); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
%legend('espaço de trabalho') 
 
6 – Determinar a trajetória do órgão efetuador (partícula A) do robô de 
2gdl, quando o motor girar de 0 a 900 graus com incremento de 1 grau e o 
fuso parte de h = 0. 
% Espaço de trabalho - 2 gdl (900 graus e deslocamento de altura) 
 
clear all 
close all 
r = 0.5; 
rd = pi/180.0; 
g = 1/rd; 
dteta = 1*rd; 
h = 0.0; 
dh = 0.0001; 
 teta = 0.0*rd; 
 for i = 1:1:900 
 z = h; 
 rav = [r*cos(teta);r*sin(teta);z]; 
 ravx(i) = rav(1); 
 ravy(i) = rav(2); 
 ravz(i) = rav(3); 
 teta = teta + dteta; 
 h =h + dh; 
end 
 plot3(ravx,ravy,ravz); grid; 
xlabel(' X (m)') 
-0.5
0
0.5
-0.5
0
0.5
0
0.2
0.4
0.6
0.8
 Z
 (
m
)
 X (m) Y (m)
ylabel(' Y (m)') 
zlabel(' Z (m)') 
%legend('trajetoria') 
 
7 – Determinar a trajetória do órgão efetuador do robô de 3 gdl, quando o 
motor girar de 0 a 900 graus com incremento de 1 grau, o fuso Z parte de 
h= 0 e o fuso Y parte de r = 0,1 m. 
 
% Espaço de trabalho - 3 gdl (900 graus e deslocamento de raio e altura) 
 
clear all 
close all 
 r = 0.1; 
dr = 0.001; 
rd = pi/180.0; 
-0.5
0
0.5
-0.5
0
0.5
0
0.02
0.04
0.06
0.08
0.1
 X (m) Y (m)
 Z
 (
m
)
g = 1/rd; 
dteta = 1*rd; 
h = 0.0; 
dh = 0.0001; 
teta = 0.0*rd; 
 
 for i = 1:1:900 
 z = h; 
 rav = [r*cos(teta);r*sin(teta);z]; 
 
 ravx(i) = rav(1); 
 ravy(i) = rav(2); 
 ravz(i) = rav(3); 
 
 teta = teta + dteta; 
 r = r+dr; 
 h =h + dh; 
end 
 
plot3(ravx,ravy,ravz,'k'); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
%legend('trajetoria') 
 
-1
-0.5
0
0.5
1
-1
-0.5
0
0.5
1
0
0.02
0.04
0.06
0.08
0.1
 X (m) Y (m)
 Z
 (
m
)
8 – Determinar o espaço de trabalho do robô esférico de dois graus de 
liberdade (2 gdl) mostrado abaixo, para e  variando de 0 a 360 graus 
com incremento de 1 grau e r = 0,5 m. 
 
% Espaço de trabalho - 2 gdl (teta e beta 360 graus) 
 
 
clear all 
close all 
 
r = 0.5; 
 
rd = pi/180.0; 
g = 1/rd; 
 
dteta = 1*rd; 
beta = 0.0*rd; 
dbeta = 1*rd; 
 
for j =1:1:360 
 
 teta = 0.0*rd; 
 
 for i = 1:1:360 
 
 rpv = [r*cos(beta)*cos(teta);r*sin(teta)*cos(beta);r*sin(beta)]; 
 
 rpvx(i,j) = rpv(1); 
 rpvy(i,j) = rpv(2); 
 rpvz(i,j) = rpv(3); 
 
 teta = teta + dteta; 
 end 
 beta = beta + dbeta; 
end 
 plot3(rpvx,rpvy,rpvz); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
 
9 – Determinar a trajetória do órgão efetuador do robô de 2 gdl, para 
e  variando de 0 a 90 graus com incremento de 1 grau e r = 0,5 m. 
% Espaço de trabalho - 2 gdl (teta e beta 90 graus) 
 
clear all 
close all 
 r = 0.5; 
rd = pi/180.0; 
g = 1/rd; 
dteta = 1*rd; 
beta = 0.0*rd; 
dbeta = 1*rd; 
 teta = 0.0*rd; 
 for i = 1:1:90 
-0.5
0
0.5
-0.5
0
0.5
-0.5
0
0.5
 Z
 (
m
)
 X (m) Y (m)
 rpv = [r*cos(beta)*cos(teta);r*sin(teta)*cos(beta);r*sin(beta)]; 
 rpvx(i) = rpv(1); 
 rpvy(i) = rpv(2); 
 rpvz(i) = rpv(3); 
 teta = teta + dteta; 
 beta = beta + dbeta; 
end 
 
plot3(rpvx,rpvy,rpvz); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
%legend('trajetoria') 
 
10 – Determinar a trajetória do órgão efetuador do robô de 3 gdl, mostrado 
abaixo, para e  variando de 0 a 60 graus com incremento de 1 grau e 
o raio inicial de 0,1 m. 
% Espaço de trabalho - 3 gdl (teta e beta 60 graus) 
 
clear all 
close all 
r = 0.1; 
dr = 0.01; 
rd = pi/180.0; 
g = 1/rd; 
0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
0.3
0.4
0
0.1
0.2
0.3
0.4
0.5
 X (m) Y (m)
 Z
 (
m
)
dteta = 1*rd; 
beta = 0.0*rd; 
dbeta = 1*rd; 
 teta = 0.0*rd; 
 for i = 1:1:60 
 rpv = [r*cos(beta)*cos(teta);r*sin(teta)*cos(beta);r*sin(beta)]; 
 rpvx(i) = rpv(1); 
 rpvy(i) = rpv(2); 
 rpvz(i) = rpv(3); 
 teta = teta + dteta; 
 r = r + dr; 
 beta = beta + dbeta; 
end 
 plot3(rpvx,rpvy,rpvz,'k'); grid; 
xlabel(' X (m)') 
ylabel(' Y (m)') 
zlabel(' Z (m)') 
%legend('trajetoria') 
 
 
 
0.1
0.15
0.2
0.25
0.3
0.35
0
0.1
0.2
0.3
0.4
0
0.2
0.4
0.6
0.8
 X (m) Y (m)
 Z
 (
m
)
ATIVIDADE 4 
 
1 – Usar o módulo simulink do Matlab, realizar a seguinte soma e 
determinar sua densidade espectral de potência. 
 Y = 10sen(10t) + 20sen(20t) + 30sen(30t) 
 
 
 
2 – Determinar a resposta de um motor cc, com função de transferência 
dada por G(s) = 5/(s2 + 2s), à uma entrada degrau, usando o simulink 
 
 
 
 
 
 
 
 
VIBRAÇÃO FORCADA AMORTECIDA DE SISTEMAS DE UM GRAU 
DE LIBERDADE (1 GDL) 
 
A Figura 1 mostra o modelo de um sistema de 1 GDL, amortecido. 
 
 
Figura 1 – Modelo de um sistema de 1 GDL amortecido 
 
Fazendo o equilíbrio dinâmico do modelo de um sistema de 1 GDL, 
amortecido mostrado na Figura 1, conforme o DCL tem-se: 
 
120 0I m a IF F F F F F u+ =  + + + + = (1) 
 
Na direção x, tem-se de (1); 
 
( ) ( ) ( ) ( )mx t cx t kx t u t+ + = (2) 
onde: 
/n k m = - frequência circular natural do sistema (rad/s); 
m – massa do sistema (kg); 
k – constante de rigidez do sistema (é função do módulo de elasticidade de 
cada sistema) (N/m); 
c – coeficiente de amortecimento viscoso (Ns/m). 
 
A transformada de Laplace de 2) é: 
 
2( ) ( ) ( )ms cs k X s U s+ + = (3) 
 
De 3), tem-se a função de transferência do sistema: 
 
2
( ) 1
( )
( ) ( )
X s
G s
U s ms cs k
= =
+ +
 (3) 
 
3 – Dada a função de transferência do sistema da Figura 1, pela eq. 3), 
determinar a resposta x(t) do sistema para u(t) = 10sen(20t). Considerando: 
m = 1 kg; c = 2 Ns/m e k = 100 N/m 
 
 
 
4 – Dada a função de transferência do sistema da Figura 1, pela eq. 3), 
determinar usando o script do Matlab a resposta x(t) do sistema para u(t) = 
10sen(20t). 
%Programa massa-mola 
clear all 
close all 
m =1; 
c = 2; 
k = 100; 
n = [1]; 
d = [1 2 100]; 
t = 0; 
dt = 0.1; 
for i = 1:1:501; 
 u = 10*sin(20*t); 
 tv(i) = t; 
 uv(i) = u; 
 t = t + dt; 
end 
y = lsim(n,d,uv,tv); 
figure(1) 
plot(tv,uv);grid 
figure(2) 
plot(tv,y);grid 
 
 
ATIVIDADE 5 
 
DETERMINAÇÃO DA VAZÃO DE UM FLUIDO USANDO EQUAÇÃO 
DE BERNOULLI 
 
O reservatório da figura abaixo é mantido com nível constante de 
fluido H, na posição 1, por uma bomba. Pede-se determinar a vazão de 
fluido no duto de diâmetro D, na posição 2. 
 
Usando a equação de Bernoulli, nos pontos 1 e 2, tem-se: 
 
2 2
1 1 2 2
1 2
2 2
 1)
v p v p
z z
g g 
+ + = + +
 
Mas, 
 
1 1 20; 01 2 = ; z ; z v p p H= = =
 Logo, 
 
2
2
2 2
2
 = 2)
v
H v gH
g
= 
 
Substituindo a velocidade obtida em 2), em 3);tem-se a vazão. 
 
2
2 2 (
4
) 3)
D
Q v A v

= =
 
1 – Elaborar um programa computacional, para obter a vazão Q do fluido 
em função da altura H. Faca H variar de 1 m até 10 m com incremento de 
1 m, para D = 0,2 m. Mostre o resultado graficamente. 
%Programa vazão 
clear all; 
close all 
 d = 0.2; 
a = (pi*d^2)/4; 
g = 9.8; 
h = 1; 
dh = 1; 
for i = 1:1:10 
 v = sqrt(2*g*h); 
 q = v*a; 
 hv(i) = h; 
 qv(i) = q; 
 h = h+dh; 
end 
figure(1) 
plot(hv,qv);grid; 
xlabel(' Altura do reservatório') 
ylabel(' Vazão') 
 
 
 
2 - Elaborar um programa computacional, para obter a vazão Q do fluido 
em função da altura H e do diâmetro D. Faca H variar de 1 m até10 m com 
incremento de 1 m e D variar de 0,2 m até 1 m com incremento de 0,1 m.. 
Mostre o resultado graficamente. 
%Programa vazão com variação do diâmetro 
clear all; close all 
 d = 0.2; deltad = 0.001; g = 9.8; h = 1; 
deltah = 1; 
 for j = 1:1:800 
 h = 1; 
for i = 1:1:10 
 v = sqrt(2*g*h); 
 a = (pi*d^2)/4; 
 q = v*a; 
 hv(i,j) = h; 
 dv(i,j) = d; 
 qv(i,j) = q; 
 h = h+deltah; 
end 
 d = d + deltad; 
end 
figure(1) 
plot(hv,qv);grid; 
1 2 3 4 5 6 7 8 9 10
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
 Altura do reservatório
 V
a
z
ã
o
xlabel(' Altura do reservatório') 
ylabel(' Vazão') 
legend('d1','d2','.....','d800') 
figure(2) 
plot3(hv,dv,qv) 
 
 
1 2 3 4 5 6 7 8 9 10
0
2
4
6
8
10
12
 Altura do reservatório
 V
a
z
ã
o
 
 
d1
d2
.....
d800
0
2
4
6
8
10
0.2
0.4
0.6
0.8
1
0
2
4
6
8
10
12
DETERMINAÇÃO DA TAXA DE CALOR EM PAREDE DE UM 
FORNO INDUSTRIAL 
 
A parede de um forno industrial é construida em tijolo refratário de 
espessura “e” e a condutividade térmica k = 1,7 W/mk. Mediçoes efetuadas 
no regime estacionário revelaram temperaturas de 1400 e 1150 k, nas 
superfícies interna e externa da parede do forno. Qual a taxa de calor 
perdida através de uma parede com dimensões de: h = 0,5 m (altura), L = 
3,0 (comprimento), e e = 0,15 m (espessura). 
 
A transmissão de calor é por condução, e a taxa de calor por unidade de 
área é dada por: 
 
int. . = -k( ) 1)extT TT
qa k
x e
−
= −

 
Substituindo os dados em 1), tem-se: 
1150 1400
)
0,15
2-1,7( 2833 W/m 2)qa
−
= = 
E, q = qa(A) = 2833(0,5)(3,0) = 4250 W 
3 - Elaborar um programa computacional no matlab, considerando T
inicial de -250 k, com incremento de 10 k e determinar a taxa de calor q. 
%Programa taxa de calor 
clear all 
close all 
h = 0.5; 
l = 3; 
e = 0.15; 
k = 1.7; 
deltat = -250; 
ddeltat = 10; 
for i = 1:1:10 
 qa = -(k*deltat)/e; 
 a = h*l; 
 q = qa*a; 
 deltatv(i) = deltat; 
 qv(i)=q; 
 deltat = deltat + ddeltat; 
end 
figure(1) 
plot(deltatv,qv);grid; 
xlabel('Variação de Temperatura') 
ylabel('Calor perdido') 
 
 
-250 -240 -230 -220 -210 -200 -190 -180 -170 -160
2600
2800
3000
3200
3400
3600
3800
4000
4200
4400
Variação de Temperatura
C
a
lo
r 
p
e
rd
id
o
4 - Elaborar um programa computacional no matlab, considerando T
inicial de -250 k, com incremento de 10 k, a espessura inicial e = 0.15 m, 
com incremento de 0.05 m, e determinar a taxa de calor q. 
%Programa taxa de calor com variação da espessura de parede 
clear all 
close all 
 
h = 0.5; 
l = 3; 
e = 0.15; 
deltae = 0.05; 
 
k = 1.7; 
ddeltat = 10; 
 
for j=1:1:5 
 
 deltat = -250; 
 
for i = 1:1:10 
 qa = -(k*deltat)/e; 
 a = h*l; 
 q = qa*a; 
 deltatv(i,j) = deltat; 
 qv(i,j)=q; 
 ev(i,j)=e; 
 deltat = deltat + ddeltat; 
end 
 e = e + deltae; 
end 
figure(1) 
plot(deltatv,qv);grid; 
xlabel('Variação de Temperatura') 
ylabel('Calor perdido') 
legend('e1', 'e2', 'e3', 'e4', 'e5') 
figure(2) 
plot3(deltatv,ev,qv);grid; 
 
 
 
 
 
 
 
 
 
 
-250 -240 -230 -220 -210 -200 -190 -180 -170 -160
1000
1500
2000
2500
3000
3500
4000
4500
Variação de Temperatura
C
a
lo
r 
p
e
rd
id
o
 
 
e1
e2
e3
e4
e5
-260
-240
-220
-200
-180
-160
0.1
0.2
0.3
0.4
0.5
1000
2000
3000
4000
5000
ATIVIDADE 6 
 
DETERMINAÇÃO DE FORÇAS EM UMA VIGA SUBMETIDA A 
CARREGAMENTO 
 
A figura A abaixo mostra uma viga submetida a um carregamento 
concentrado. Pede-se determinar a reação no rolete B. 
Do Diagrama de Corpo Livre (DCL) da Figura B, tem-se: 
 
0
/aP L
 =  
= −
 A B B
B
 M =0 AC x P + AB x R a i x -P j+L i x R j = 0
 R
 
 
 
 
 
1 - Elaborar um programa computacional no matlab, considerando P = 
1000 N, L = 0.8 m, “a” inicial igual a 0,2 m, com incremento de 0,1 m e 
determinar a reação no apoio B. 
 
%Programa viga 
 
clear all 
close all 
 
a = 0.2; 
deltaa = 0.1; 
p = 1000.0; 
l = 0.8; 
 
for i = 1:1:7; 
 av = [a;0;0]; 
 pv = [0;-p;0]; 
 m1 = cross(av,pv); 
 rby = -m1(3)/l; 
 rbv = [0;rby;0]; 
 rav = -rbv-pv; 
 avv(i) = av(1); 
 rbvv(i) = rbv(2); 
 a = a + deltaa; 
end 
figure(1) 
plot(avv,rbvv);grid 
xlabel(' Variação do ponto de aplicação da carga') 
ylabel(' Reação no apoio B') 
 
 
 
2 – Elaborar um programa computacional no matlab, considerando P inicial 
igual a 1000 N, com incremento de 100 N, L = 0.8 m, “a” inicial igual a 
0,2 m, com incremento de 0,1 m e determinar a reação no apoio B. 
 
%Programa viga com variação da carga 
 
clear all 
close all 
 
a = 0.2; 
deltaa = 0.1; 
p = 1000.0; 
deltap = 100.0; 
l = 0.8; 
for j =1:1:5 
 a = 0.2; 
for i = 1:1:7; 
 av = [a;0;0]; 
 pv = [0;-p;0]; 
 m1 = cross(av,pv); 
 rby = -m1(3)/l; 
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
200
300
400
500
600
700
800
900
1000
 Variação do ponto de aplicação da carga
 R
e
a
ç
ã
o
 n
o
 a
p
o
io
 B
 rbv = [0;rby;0]; 
 rav = -rbv-pv; 
 avv(i,j) = av(1); 
 rbvv(i,j) = rbv(2); 
 a = a + deltaa; 
end 
p = p + deltap; 
end 
figure(1) 
plot(avv,rbvv); grid 
xlabel(' Variação do ponto de aplicação da carga') 
ylabel(' Reação no apoio B') 
legend('P1','P2','P3','P4','P5') 
 
 
 
 
 
 
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
200
400
600
800
1000
1200
1400
 Variação do ponto de aplicação da carga
 R
e
a
ç
ã
o
 n
o
 a
p
o
io
 B
 
 
P1
P2
P3
P4
P5
A figura abaixo mostra uma viga engastada na extremidade direita e 
livre na esquerda. Na extremidade livre atua uma carga concentrada W. A 
equação da deflexão “y” da viga é dada abaixo. 
 
 ( )− +3 2 3W
 y = x 3L x-2L
6EI
 
 
3 – Elaborar um programa computacional, para obter a deflexão “y” da 
viga em função do comprimento x . Faca x variar de 0 m até 1 m com 
incremento de 0,1 m, considerando: E = 2.109 N/m2 , I = 0,5 m4 , W = 3000 
N. Mostre o resultado graficamente. 
 
%Programa deflexão da viga 
 
l = 1.0; %em m 
me = 2e9; %módulo de elasticidade 
mi = 0.5; %momento de inércia 
w = 3000.0; %carga concentrada 
deltax = 0.1; 
 
 x = 0.0; 
for j = 1:1:11; 
 y =(w/(6*me*mi))*(-x^3+3*l^2*x-2*l^3); 
 xv(j) = x; 
 yv(j) = y; 
 miv(j) = mi; 
 x = x+deltax; 
end 
 
 
figure(1) 
plot(xv,yv);grid; 
xlabel('comprimento (m)') 
ylabel('deflexão (m)') 
 
 
 
4 – Elaborar um programa computacional, para obter a deflexão “y” da 
viga em função do comprimento x e do momento de inércia I. Faca x variar 
de 0 m até 1 m com incremento de 0,1 m, I variar de 0,5 m4 até 0,8 m4, 
considerando: E = 2.109 N/m2 , I = 0,5 m4 , W = 3000 N. Mostre o 
resultado graficamente. 
 
%Programa deflexão da viga com variação do momento de inércia 
 
clear all 
close all 
 
l = 1.0; %em m 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0
x 10
-6
comprimento (m)
d
e
fl
e
x
ã
o
 (
m
)
me = 2e9; %módulo de elasticidade 
mi = 0.5; %momento de inércia 
w = 3000.0; %carga concentrada 
deltax = 0.1; 
deltami = -0.1; 
 
for i = 1:1:3; 
 x = 0.0; 
for j = 1:1:11; 
 y =(w/(6*me*mi))*(-x^3+3*l^2*x-2*l^3); 
 xv(j,i) = x; 
 yv(j,i) = y; 
 miv(j,i) = mi; 
 x = x+deltax; 
end 
mi = mi+deltami; 
end 
 
figure(1) 
plot(xv,yv);grid; 
xlabel('comprimento(m)') 
ylabel('deflexão(m)') 
 
figure(2) 
plot3(xv,miv,yv);grid; 
xlabel('comprimento(m)') 
ylabel('momento de inércia(m^4)') 
zlabel('deflexão(m)') 
 
 
 
 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1.8
-1.6
-1.4
-1.2
-1
-0.8
-0.6
-0.4
-0.2
0
x 10
-6
comprimento(m)
d
e
fl
e
x
ã
o
(m
)
0
0.2
0.4
0.6
0.8
1
0.35
0.4
0.45
0.5
-2
-1.5
-1
-0.5
0
x 10
-6
comprimento(m)momento de inércia(m4)
d
e
fl
e
x
ã
o
(m
)
A figura abaixo mostra uma viga engastadana extremidade direita e 
livre na esquerda, submetida a um carregamento distribuído w. A equação 
da deflexão “y” da viga é dada abaixo. 
 
 ( )− +
4 3 4w x L x L
 y = - 
E I 24 6 8
 
 
 
5 – Elaborar um programa computacional, para obter a deflexão “y” da 
viga em função do comprimento x. Faca x variar de 0 m até 1 m com 
incremento de 0,1 m, considerando: E = 200.109 N/m2 , I = 0,5 m4 , W = 
3000 N. Mostre o resultado graficamente. 
 
 
%Programa viga carregamento distribuído 
 
clear all 
close all 
 
l = 1.0; %em m 
me = 200e9; %módulo de elasticidade 
mi = 0.5; %momento de inércia 
w = 3000.0; %carga concentrada 
deltax = 0.1; 
 x = 0.0; 
for j = 1:1:11; 
 y =(w/(me*mi))*((-x^4/24)+((l^3)*x/6)-(l^4/8)); 
 xv(j) = x; 
 yv(j) = y; 
 miv(j) = mi; 
 x = x+deltax; 
end 
 
 
figure(1) 
plot(xv,yv);grid; 
xlabel('comprimento (m)') 
ylabel('deflexão (m)') 
 
 
 
 
6 - Elaborar um programa computacional, para obter a deflexão “y” da viga 
em função do comprimento x e de w. Faca x variar de 0 m até 1 m com 
incremento de 0,1 m, sendo w inicial igual a 3000 N e deltaw = 1000 N, 
considerando: E = 200.109 N/m2 , I = 0,5 m4 . Mostre o resultado 
graficamente. 
 
 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
x 10
-9
comprimento (m)
d
e
fl
e
x
ã
o
 (
m
)
%Programa viga carregamento distribuído e variação da carga 
 
clear all 
close all 
 
l = 1.0; %em m 
me = 200e9; %módulo de elasticidade 
mi = 0.5; %momento de inércia 
w = 3000.0; %carga distribuida 
deltaw = 1000; 
deltax = 0.1; 
 
for j = 1:1:6 
 x = 0.0; 
for i = 1:1:11; 
 y =(w/(me*mi))*((-x^4/24)+((l^3)*x/6)-(l^4/8)); 
 xv(i,j) = x; 
 yv(i,j) = y; 
 wv(i,j) = w; 
 x = x+deltax; 
end 
w = w + deltaw; 
end 
 
figure(1) 
plot(xv,yv);grid; 
xlabel('comprimento (m)') 
ylabel('deflexão (m)') 
 
figure(2) 
plot3(xv,wv,yv);grid; 
 
 
 
 
 
 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0
x 10
-8
comprimento (m)
d
e
fl
e
x
ã
o
 (
m
)
0
0.2
0.4
0.6
0.8
1
2000
4000
6000
8000
-1
-0.8
-0.6
-0.4
-0.2
0
x 10
-8
ATIVIDADE 7 
 
 
O movimento retilineo uniformemente variado (mruv) é descrito 
pelas equações abaixo: 
 
 
2
0 0 0,5
at
x x v t at
+
= + +
0
 a = k ( constante)
v = v 
 
 
1 - Faça a = 10 m/s2 , v0 = 0 e x0 = 0 e considere t inicial igual a zero, com 
incremento de 0,1 seg e contrua os gráficos de a, v e x em função do tempo 
t. Grave o tempo, a velocidade e o deslocamento em um arquivo de dados. 
 
 
%Programa mruv 
 
clear all 
close all 
 
 
a = 1; 
vo = 0; 
xo = 0; 
 
t = 0; 
dt = 0.1; 
 
 
for i = 1:1:50; 
 a = 10; 
 v = vo+a*t; 
 x = xo+vo*t+0.5*a*t^2; 
 
 xv(i) = x; 
 av(i) = a; 
 vv(i) = v; 
 tv(i) = t; 
 t = t + dt; 
end 
 
figure(1) 
plot(tv,av,tv,vv,tv,xv);grid; 
legend('a','v','x') 
figure(2) 
plot3(tv,vv,xv);grid; 
 
 resul = [tv' vv' xv']; 
 
save resultado.dat resul -ascii 
 
 
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
20
40
60
80
100
120
140
 
 
a
v
x
 
 
2 – Faça a leitura do arquivo de dados do programa anterior e construa o 
gráfico da velocidade e do deslocamento em função do tempo. 
 
% Leitura de arquivo de dados 
 
clear all 
close all 
 
load resultado.dat 
 
 t = resultado(:,1); 
 
 v = resultado(:,2); 
 
 x = resultado(:,3); 
 
figure(1) 
plot(t,v,t,x);grid; 
 
figure(2) 
plot3(t,v,x);grid; 
 
0
1
2
3
4
5
0
20
40
60
0
50
100
150
GEOMETRIA DO MOVIMENTO DO MECANISMO CURSOR-
MANIVELA 
 
A Figura 2.1, mostra as linhas de centro das peças do mecanismo 
cursor-manivela. O objetivo é a determinação da posição angular da biela 
(peça 3) e a posição linear do centro do cursor (peça 4), em função da 
posição angular da manivela (peça 2), e dos comprimentos O2A = r2, da 
manivela e AB = r3, da biela. 
 
A Figura 2.1 – Mecanismo cursor - manivela 
 
 A Figura 2.2, mostra o polígono de vetores construído a partir da 
Fig. 2.1, que será usado para obtenção da posição angular da biela (peça 3) 
e a posição linear do centro do cursor (peça 4). 
 
Figura 2.2 – Polígono de vetores do mecanismo cursor - manivela 
 
 A solução é obtida vetorialmente, considerando os vetores e a base 
ortonormal da Fig. 2.2. Da adição dos vetores, tem-se: 
 
 r2 + r3 = rB 2) 
 
Escrevendo 2) em função das componentes de cada vetor, tem-se: 
 
 










=










+










0
0
0
)(
)cos(
0
)(
)cos(
33
33
22
22 Br
senr
r
senr
r
 




 3) 
 
Separando as componentes de 3) segundo os eixos X e Y, tem-se: 
 
 
4) BrrrX =+ )cos()cos(: 3322 
 
 
 0)()(: 3322 =+  senrsenrY 5) 
 
Do sistema de equações 4) e 5), tem-se da Eq. 6): 
 
 ))( 2
3
2
3  sen
r
r
(-sen 1-= 6) 
 
Substituindo ϴ3 obtido em 6) em 4), tem-se: 
 
 7) )cos()cos( 3322  rrrB += 
3 – Elaborar um programa computacional para análise da geometria do 
movimento do mecanismo cursor – manivela, em um ciclo de movimento 
da manivela (peça 2). Mostrar os resultados graficamente. 
% Programa Cursor- Manivela 
 
clear all; 
close all; 
format long e 
 r2 = 0.10; 
r3 = 0.20; 
 convrd=pi/180; 
convg=1/convrd; 
teta2=0.0*convrd; 
dteta2 = 0.5*convrd; 
 
 
 for i = 1:1:720; 
 r2v = [r2*cos(teta2);r2*sin(teta2);0]; 
 teta3 = asin((-r2v(2)/r3)); 
 r3v = [r3*cos(teta3);r3*sin(teta3);0]; 
 rbv = r2v + r3v; 
 teta3v(i) = teta3*convg; 
 rbvx(i) = rbv(1); 
 teta2v(i) = teta2*convg; 
 teta2 = teta2 + dteta2; 
 end 
 
 % curso do cursor 
 
 sx = max(rbvx)-min(rbvx) 
 
figure(1) 
plot(teta2v,teta3v);grid 
xlabel('teta2 (graus)') 
ylabel('teta3 (graus)') 
 
figure(2) 
plot(teta2v,rbvx);grid 
xlabel('teta2 (graus)') 
ylabel('rb (m)') 
 
 
 
0 50 100 150 200 250 300 350 400
-40
-30
-20
-10
0
10
20
30
teta2 (graus)
te
ta
3
 (
g
ra
u
s
)
0 50 100 150 200 250 300 350 400
0.1
0.15
0.2
0.25
0.3
0.35
teta2 (graus)
rb
 (
m
)
4 – Elaborar um programa computacional para análise da geometria do 
movimento do mecanismo cursor – manivela, em um ciclo de movimento 
da manivela (peça 2). Faça também a variação no comprimento da 
manivela. Mostrar os resultados graficamente. 
 
% Programa Cursor- Manivela com variação de comprimento da manivela 
 
clear all; 
close all; 
format long e 
 r2 = 0.10; 
 dr2 = 0.01; 
r3 = 0.20; 
 convrd=pi/180; 
convg=1/convrd; 
teta2=0.0*convrd; 
dteta2 = 0.5*convrd; 
 
 for j = 1:1:6 
 
teta2=0.0*convrd; 
 for i = 1:1:720; 
 r2v = [r2*cos(teta2);r2*sin(teta2);0]; 
 teta3 = asin((-r2v(2)/r3)); 
 r3v = [r3*cos(teta3);r3*sin(teta3);0]; 
 rbv = r2v + r3v; 
 teta3v(i,j) = teta3*convg; 
 rbvx(i,j) = rbv(1); 
 teta2v(i,j) = teta2*convg; 
 teta2 = teta2 + dteta2; 
 end 
 r2 = r2 + dr2 
 end 
 % curso do cursor 
 
 sx = max(rbvx)-min(rbvx) 
 
figure(1) 
plot(teta2v,teta3v);grid 
xlabel('teta2 (graus)') 
ylabel('teta3 (graus)') 
 
figure(2) 
plot(teta2v,rbvx);grid 
xlabel('teta2 (graus)') 
ylabel('rb (m)') 
 
 
 
 
 
0 50 100 150 200 250 300 350 400
-50
-40
-30
-20
-10
0
10
20
30
40
50
teta2 (graus)
te
ta
3
 (
g
ra
u
s
)
0 50 100 150 200 250 300 350 400
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
teta2 (graus)
rb
 (
m
)
ATIVIDADE 8 
 
REPRESENTAÇÃO DE SISTEMAS DINÂMICOS NO ESPAÇO DE 
ESTADOSUm sistema dinâmico pode ser representado por um modelo, através 
de equações diferenciais ordinárias, onde o tempo é a variável 
independente. 
Uma equação diferencial de ordem n pode ser representada por uma 
equação matricial-vetorial diferencial de primeira ordem. Se os n elementos 
do vetor são um conjunto de variáveis de estado, então a equação matricial-
vetorial diferencial é chamada de equação de estado. 
 
REPRESENTAÇÃO NO ESPAÇO DE ESTADOS DE SISTEMAS 
DESCRITOS POR EQUAÇÕES DIFERENCIAIS LINEARES DE 
ORDEM n SEM DERIVADAS DE EXCITAÇÃO 
 
Considere o sistema representado pelo modelo de ordem n abaixo: 
 
 
 1) 
 
 
De 1); 
 
 
2) 
 
 
)()(
)(
....
)()(
011
1
1 tgtya
dt
tdy
a
dt
tyd
a
dt
tyd
n
n
nn
n
=++++
−
−
−
)()(
)(
...
)()(
011
1
1 tgtya
dt
tdy
a
dt
tyd
a
dt
tyd
n
n
nn
n
+−−−−=
−
−
−
Definindo n (número da ordem da equação diferencial original) 
novas variáveis, x1(t), x2(t),...., xn(t), pelas equações: 
 
1
1
2
2
321
)(
)(,
)(
)(,
)(
)(),()(
−
−
====
n
n
n
dt
tyd
tx
dt
tdy
tx
dt
tdy
txtytx ,.......... 3) 
 
As novas variáveis de 3) são interligadas pelas equações 4). 
 
 










=
=
=
=
− )()(
..................
..................
)()(
)()(
)()(
1
43
32
21
txtx
txtx
txtx
txtx
nn




 4) 
 
Escrevendo a Eq. 2) em função das novas variáveis dadas por 3); 
 
 )()()(...)(
)(
10211 tgtxatxatxa
dt
tyd
nnn
n
+−−−−= −
 5) 
 
Diferenciando a última equação de 3); 
 
n
n
n
n
nn
n
n
dt
tyd
dt
tyd
dt
d
tx
dt
tyd
tx
)(
]
)(
[)(
)(
)(
1
1
1
1
===
−
−
−
−
  6) 
 
Usando 6) em 5); 
 
 )()()(...)( 10211 tgtxatxatxax nnn +−−−−= −
 7) 
 
Adicionando 7) em 4); 
 
 











+−−−−=
=
=
=
=
−
−
)()()(...)(
)()(
..................
..................
)()(
)()(
)()(
10211
1
43
32
21
tgtxatxatxax
txtx
txtx
txtx
txtx
nnn
nn





 8) 
 
 
Escrevendo 8) na forma matricial: 
 
 
)(
1
0
.
.
.
0
0
0
.
.
.
...
10..0000
........
........
........
00..1000
00..0100
00..0010
.
.
.
1
3
2
1
13210
1
3
2
1
tg
x
x
x
x
x
aaaaax
x
x
x
x
n
n
nn
n


























+




















































−−−−
=


























−
−
−





 9) 
 
 
 Ou de 9): 
 
 
( ) ou,X AX Bg t
X AX Bu
= +
= +
 10) 
 
Onde: X – Vetor de Estados 
 A – Matriz Dinâmica 
 B – Matriz de Entrada 
 u – entrada ou excitação 
 
A saída ou resposta Y(t) (variável original) é dada por: 
 
 Y(t) = Cx + Du 11) 
 
C e D – Matrizes de Saída 
 
 
VIBRAÇÃO LIVRE E FORCADA AMORTECIDA DE SISTEMAS DE 
UM GRAU DE LIBERDADE (1 GDL) 
 
A Figura 1 mostra o modelo de um sistema de 1 GDL, amortecido. 
 
 
Figura 1 – Modelo de um sistema de 1 GDL amortecido 
Fazendo o equilíbrio dinâmico do modelo de um sistema de 1 GDL, 
amortecido mostrado na Figura 1, conforme o DCL tem-se: 
 
120 0I m a IF F F F F F u+ =  + + + + = (1) 
 
Na direção x, tem-se de (1); 
( ) ( ) ( ) ( )
( ( ) ( ) ( ))
( )
mx t cx t kx t u t
kx t cx t u t
x t
m
+ + = 
− − +
=
 (2) 
onde: 
/n k m = - frequência circular natural do sistema (rad/s); 
m – massa do sistema (kg); 
k – constante de rigidez do sistema (é função do módulo de elasticidade de 
cada sistema) (N/m); 
c – coeficiente de amortecimento viscoso (Ns/m). 
 
Transformar a Eq. 2) para a forma de espaço de estados e mostrar as 
matrizes A, B, C e D. 
 
 A equação de estados é da forma: 
 X AX Bu= + 
Definindo as variáveis de estado como segue: 
 
1 2 e x x x x= = (3) 
 
Usando (3) e (2), tem-se: 
 
 
1 2
2 1 2
1
x x
k c
x x x u
m m m
=


= − − +

 (4) 
 
Escrevendo 4) na forma matricial, tem-se: 
 
 
 
 
Onde: 
 
0 1 0
; 1
 
 
 
A Bk c
m m m
   
   = =
   − −
      
 
 
A saída é dada por: 
 
Y CX Du= + 
 
Ou seja; 
 
   1 1
2 2
1 0 1 0
0 0
0 1 0 1
 
 
y x x
u u
y x x
        
= + = +        
        
 
 
Onde: 
 
1 0
[0]
0 1
 
 e 
 
C D
 
= = 
 
 
 
 
 
1 - Dadas as matrizes A, B, C e D, do sistema da Figura 1, do item 2, 
elaborar um programa computacional, e determinar: 
 
a) A resposta forçada x(t) do sistema para u(t) = 10sen(20t); 
 
b) A resposta forçada x(t) do sistema para u(t) = 10sen(10t); 
 
 1 1
2 2
0 1 0
1
 
 
x x
uk c
x x
m m m
   
      = +      − −         
c) A resposta livre x(t) do sistema com amortecimento; 
 
d) A resposta livre x(t) do sistema sem amortecimento; 
 
Considerando: m = 1 kg; c = 2 Ns/m e k = 100 N/m 
 
%Programa massa-mola-amortecedor (letra a) 1gdl 
clear all 
close all 
 
m = 1.0; 
ca = 2.0; 
k = 100.0; 
 
a = [0 1;-k/m -ca/m]; 
 
b = [0;1/m]; 
 
c = [1 0]; 
 
d = [0]; 
 
xo = [0.1 0]; 
 t = 0.0; 
dt = 0.05; 
 for i = 1:1:500; 
 u = 10*sin(20*t); 
 tv(i) = t; 
 uv(i) = u; 
 t = t+dt; 
end 
 
y = lsim(a,b,c,d,uv,tv,xo); 
 
figure(1) 
plot(tv,uv);grid; 
figure(2) 
plot(tv,y);grid; 
 
ATIVIDADE 9 
 
A Figura 1 mostra um Sistema de 2 GDL com uma força de 
excitação nas duas massas. 
 
 
Figura 1 – Sistema de 2 GDL 
 
 
Do DCL do sistema, tem-se: 
 
1 1 1 2 1
2 2 1 2 2
2
2
mx cx kx kx u
mx cx kx kx u
+ + − =


 + − + =
 (1) 
 
Transformar (1) para a forma de espaço de estados e mostrar as matrizes A, 
B, C e D. 
Representando (1) na forma de Espaço de Estados, com a mudança de 
variáveis como (2); 
 
1 1 2 1 3 2 4 2y x y x y x e y x= = = = , , (2) 
 
Tem-se de (1): 
 
1 2
2 1 2 3 1
 
2 1
y y
k c k
y y y y u
m m m m
=
= − − + +
 
 
3 4
4 1 3 4 2
2 1
y y
k k c
y y y y u
m m m m
=
= − − +
 (3) 
 
Escrevendo 3) na forma matricial, chega-se em 4). 
 
1 1
2 2 1
3 3 2
4 4
0 1 0 0
0 02
0
1/ 0
0 0 0 1 0 0
2 0 1/
0
y yk c k
y y umm m m
y y u
k k c my y
m m m
 
      − −
      
       = +                
− −         
  
 4) 
 
A saída é dada por: Z = CY, dada por 5). 
 
1 1
2 2
3 3
4 4
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
z y
z y
z y
z y
    
    
    =
    
    
       
 5) 
 
 
1 - Dadas as matrizes A, B, C e D, do sistema da Figura 1, do item 2, 
elaborar um programa computacional, e determinar: 
 
a) A resposta forçada x(t) do sistema para u1(t) = 10sen(20t) e 
u2(t) = 5sen(15t); 
 
b) A resposta livre x(t) do sistema com amortecimento; 
 
c) A resposta livre x(t) do sistema sem amortecimento; 
 
Considerando: m = 1 kg; c = 2 Ns/m e k = 100 N/m 
 
%Programa massa-mola-amortecedor (letra a) 2gdl 
 
clear all 
close all 
 
m = 1.0; 
ca = 2.0; 
k = 100.0; 
 
a = [0 1 0 0;(-2*k)/m -ca/m k/m 0;0 0 0 1;k/m 0 (-2*k)/m -ca/m]; 
b = [0 0;1/m 0;0 0;0 1/m]; 
c = [1 0 0 0;0 0 1 0]; 
 
c1 = [1 0 0 0]; 
c2 = [0 0 1 0]; 
d = [0];xo = [0.1 0 0.2 0]; 
 
t = 0.0; 
dt = 0.05; 
 
for i = 1:1:500; 
 u1 = 10*sin(20*t); 
 u2 = 5*sin(15*t); 
 tv(i) = t; 
 u1v(i) = u1; 
 u2v(i) = u2; 
 t = t+dt; 
end 
 u = [u1v;u2v]; 
 
 y = lsim(a,b,c,d,u,tv,xo); 
y1 = lsim(a,b,c1,d,u,tv,xo); 
y2 = lsim(a,b,c2,d,u,tv,xo); 
 
figure(1) 
plot(tv,u1v,tv,u2v);grid; 
 
figure(2) 
plot(tv,y);grid; 
legend('x1', 'x2') 
 
figure(3) 
plot(tv,y1,tv,y2);grid; 
 
 
 
 
ATIVIDADE 10 
 
A Figura 1 mostra um pêndulo simples e seu DCL. Analisando as 
forças na direção tangencial, tem-se: 
 
0IF F+ =  
 
: 0 0It Wsen F Wsen mL  − − =  − − =  
 
0mL Wsen + =  
 
 
 0
g
sen
L
 + = (1) 
 
 
 
Figura 1 - Pendulo simples 
 
A equação (1) é uma EDM de 2a ordem, de coeficientes constantes, 
homogênea e não linear. Linearizando (1); ou seja, considerando pequenos 
deslocamentos angulares, tem-se: 
 
sen  (2) 
 
Substituindo (2) em (1); 
 
0
g
L
 + = (3) 
 
E, 
2
2n
n
g L
e
L g

   

= =  = (4) 
 
A EDM (3) é de 2a ordem, de coeficientes constantes, homogênea e linear. 
Transformar a Eq. 1) e a Eq. 3) para a forma de espaço de estados. 
Definindo como variáveis de estados: 
 
(1) e y(2) 5)y  = =
 
 
Usando 5) e 1), tem-se: 
 
 
(1) (2)
 
(2) = - ( (1))
 
y y
g
y sen y
L

 =





 6) 
Usando 5) e 3), tem-se: 
 
(1) (2)
 
(2) = - ( (1))
 
y y
g
y y
L

 =





 7) 
 
1 - Elaborar um programa computacional para resolver 6) e 7), 
considerando g/L = 20,0. 
 
%Função 'se' das derivadas (linear) 
function dx=se(t,y) 
dx=[y(2);-20*y(1)]; % g/L = 20 
 
%Função 'se1' das derivadas (não linear) 
function dx=se1(t,y) 
dx=[y(2);-20*sin(y(1))]; 
 
%Programa pêndulo 
clear all 
close all 
 t0 = 0; 
dt = 0.1; 
tf = 2; 
 x0 = [0.8 0]; 
 [t y1] = ode45('se', [t0:dt:tf], x0); % modelo linear 
 [t y2] = ode45('se1', [t0:dt:tf], x0); % modelo não linear 
 figure(1) 
plot(t,y1(:,1),'k',t,y2(:,1),'r') 
legend('Solução Linear','Solução Não Linear') 
 
 
 
 
 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
 
 
Solução Linear
Solução Não Linear
ATIVIDADE 11 
 
Transformar a Eq. (1) para a forma de espaço de estados e mostrar as 
matrizes A, B, C e D. 
 100 20x x x u+ + = 
(1) 
 
Representando (1) na forma de Espaço de Estados, com a mudança 
de variáveis como (2); 
 
1 2 3 4x x x x x x e x x= = = = , , (2) 
 
Tem-se de (1): 
 
1 2
2 3
3 4
4 1 3
 
20 100
x x
x x
x x
x x x u
=

=

=
 = − − +
 (3) 
 
Escrevendo 3) na forma matricial, chega-se em 4). 
 
 
1 1
2 2
3 3
4 4
0 1 0 0 0 
0 0 1 0 0
0 0 0 1 0 
20 0 100 0 1
x x
x x
u
x x
x x
      
      
      = +
      
      
− −         
 4) 
 
A saída é dada por: Y = CX, dada por 5). 
 
1 1
2 2
3 3
4 4
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
y x
y x
y x
y x
    
    
    =
    
    
       
 5) 
 
1 - Dadas as matrizes A, B, C e D, das equações 4) e 5), do item 1, elaborar 
um programa computacional, considerando nulas as condições iniciais, 
determinar: 
 
a) A resposta forçada x(t) de 1) para u(t) = 10sen(2t); 
b) A resposta livre x(t) de 1). 
 
%Programa EDO de quarta ordem 
 
clear all 
close all 
 
a = [0 1 0 0;0 0 1 0;0 0 0 1;-20 0 -100 0]; 
 b = [0;0;0;1]; 
 c = [1 0 0 0]; 
d = [0]; 
 
xo = [0 0 0 0]; 
 
t = 0.0; 
dt = 0.1; 
 
for i = 1:1:140; 
 u = 10*sin(2*t); 
 
 tv(i) = t; 
 uv(i) = u; 
 t = t+dt; 
end 
 
 y = lsim(a,b,c,d,uv,tv,xo); 
 figure(1) 
plot(tv,uv);grid; 
 figure(2) 
plot(tv,y);grid; 
legend('x') 
Resposta a) 
 
 
 
0 2 4 6 8 10 12 14
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
 
 
x
0 5 10 15 20 25 30
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
 
 
x