Universidade Federal Fluminense

Escola de Engenharia Industrial Metalúrgica de Volta Redonda

Pólo Universitário de Volta Redonda

VIBRAÇÕES

Solução da Equação de Movimento pelo Método de Runge -

Kutta de 4ª Ordem (RK- 4)

Anderson Zenken Nakazato

Janaan Pereira da Cunha
Renner Egalon Pereira

Volta Redonda, 13 de Janeiro de 2014

Introdução e Objetivos

Este trabalho visa à obtenção da solução das equações de movimento de

um sistema massa-mola-amortecedor com 2 graus de liberdade (2 GDL’s)
numericamente, aplicando o Método de Runge – Kutta de 4ª Ordem (RK-4).

Uma vez determinadas as variáveis de estado do sistema (deslocamentos

e velocidades) através da aplicação do método numérico RK – 4, uma análise
gráfica (como forma de interpretação dos resultados) será feita envolvendo as
variáveis de estado ao longo do tempo a fim de estudar como se comportam e
interagem entre si os dois corpos do sistema, como por exemplo, estabelecer o
que seriam condições de ressonância e batimento de um sistema.

As etapas da análise dinâmica do sistema são: modelagem física,

modelagem matemática, solução do modelo matemático através do método
numérico RK – 4 e interpretação dos resultados.

A modelagem física visa representar todas as características

importantes do sistema, com o objetivo de apresentar as equações que
governam o movimento do mesmo.

Na modelagem matemática é feita a dedução das equações diferenciais

que governam o sistema mecânico, utilizando para isso, o princípio físico das
Equações de Lagrange.

A etapa de solução do modelo matemático consiste na aplicação do

método Runge – Kutta de 4ª Ordem para resolver o sistema de equações
diferenciais. As equações do sistema (de 2ª ordem) estão acopladas entre si,
com as variáveis dependentes e suas derivadas aparecem em mais de uma
equação.

A interpretação dos resultados visa comparar os resultados já esperados

teoricamente com os resultados do modelo numérico, analisando como os
corpos interagem entre si.
Fundamentos Teóricos
Modelagem Física

O sistema massa-mola-amortecedor é apresentado a seguir:

Figura 1 – Modelo Físico

Um disco de massa “ ” raio “ ” tem seu centro (ponto O) conectado a

uma mola linear que está presa a uma plataforma móvel de massa “ ”
submetido ao um forçamento harmônico dado por = cos ( ) que, por sua
vez, está conectada a outra mola linear que está presa a um referencial
inercial, conforme mostrado na figura 1. O disco rola sem deslizar sobre a
plataforma móvel e a plataforma móvel desliza sem atrito sobre a superfície
horizontal.

De acordo com a figura 2 abaixo, verifica-se que existe uma

compatibilidade geométrica entre os deslocamentos linear e angular do disco:

Figura 2 - Compatibilidade geométrica entre disco e plataforma

 2 − 1 = 

 =
2 − 1

 =
2 − 1

Modelagem Matemática

Utilizou-se o princípio físico das Equações de Lagrange para determinar

as equações de movimento. Logo abaixo são apresentadas as equações do
sistema:

− + = ( )

( �çã � � )
Onde: = � − � ∶ � � � ( . � é � �� � − . �� � � )

=
1
; = 1
; cos 
=
2 2 0
Energias do sistema:
Plataforma:
1 2
� = 1 ( �� � é � � � � �çã )
2
Disco:
1 2
� = 2 ( �� � é � � � � �çã )
2
1 1 ² 2 − 1 ² 1
0� −
2
= = = 2 1 ²
2 2 2 ² 4

( �� � é � � �çã )
Molas:
 1
= ( 2 − 1 )²
2
( �� �� á � �)
Amortecedor:
1
= 2 − 1 ²
2
( �� � � � �)

Após a obtenção das energias do sistema e posteriormente utilizadas na

Equação de Lagrange, tem-se as equações de movimento do sistema:
1 2 1 2 1 1 1
= 1 + 2 + 2 − 1 ²− 2 − 1
2
− 2 − 1 ²
2 2 4 2 2
( � � � )

+ 1 − 2 + 2 1 − 2 + 1 − 2 = cos ( )
2 2

3
2 − 1 − 1 + 2 − 1 + 2 =0
2 2

Na forma matricial:

( )
+ −
2 2 1 − 1 2 − 1 cos⁡
+ + =
3 − − 2 0
− 2 2
2 2
Solução do Modelo Matemático pelo Método RK- 4

O método de Runge-Kutta só resolve sistema de equações diferenciais

de primeira ordem, portanto é necessário transformar as equações de
movimento, dadas naturalmente por equações diferenciais de segunda ordem
em um sistema de equações de primeira ordem, explicitando as variáveis de
estado, conforme a seguir. O sistema de equações discretizadas, segundo o
método RK-4 para a resolução numérica também é mostrado abaixo.

1 =
1 = =
1 = =
Variáveis de estado
2 =
2 = =
2 = =

2 −5 + 2 − +3 
3 +
=

2 − −2 + 2 − + ( )
(3 + )

Método de Runge-Kutta 4ª Ordem:

∆
� +1
= � + ( 1 +2 2 +2 3 + 4)
6

1 = � ,

∆ ∆
2 = � + 1 , +
2 2
∆ ∆
3 = � + 2 , +
2 2
 4 = � + 3 ∆ , +∆

+1
+1
� +1
= +1
� = =
+1

11 21 31 41
12 22 32 42
1 = 2 = 3 = 4 =
13 23 33 43
14 24 34 44

O código computacional do Runge-Kutta 4ª Ordem foi programado

utilizando o MATLAB, um software interativo para cálculo numérico. O código é
dividido em três “scripts”, o primeiro com o código do RK-4 efetivamente, o
segundo com o programa principal e o terceiro onde se entram com as funções
(equações diferenciais de 1ª ordem discretizadas). O código se encontra
anexado ao final do relatório.

Resultados e Discussão

Adotaram-se os seguintes parâmetros para o estudo inicial do sistema

mecânico:

Tabela 1 – Parâmetros do sistema

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 1 N.s/m

Frequência de forçamento  1 rad/s

Amplitude de forçamento � 1N
 Para calcular as frequências e os modos de vibração, é possível fazer
uma análise de vibração livre, calculando as frequências naturais e modos de
vibração considerando a matriz de amortecimento e forçamento nulas. A
solução recai num problema de autovalores e autovetores, conforme a seguir:

−  = 0

−
+ − 1 0
2 − 2 2 =
− 3
− 2 0
2 2

Fazendo  = ²

2 −  + − + 
3
2 2
− + 
=0
−
2 2

2 +
3 2
−5
+ − + 2
=0
2 2 2

(Polinômio Característico)

5 17
1,2 =
+ ± ²− + ²
2 4
3 + ²

Utilizando os parâmetros adotados para o sistema (tabela 1) e

substituindo na expresão acima, tem-se:

1 =
1
5

2 = 1

Logo, as frequências naturais são:

1 =
5
� / = 0,447213595 � /
5
2 = 1 � /
 Retornando à expressão original e substituindo novamente os
parâmetros adotados, calculam-se os autovalores (modos de vibração):

2 −  + − +  0
3
1
2 2
− + 
=
− 2 0
2 2

� �= :
1
5

8 4
− 1 0
5 5 =
4 2
− − 2 0
5 5

8 4
1 − 2 =0 2 =2 1
5 5

4 2
− 1 + 2 =0 2 =2 1
5 5

1
1 =
2

Percebe-se que no primeiro modo 1, quando a plataforma se desloca

de 1 o disco desloca o dobro, de 2.

� �  = 1:
1 0
0 0
=
0 −2
2 0

0+0 =0
−2 2 =0 2 =0
 2 = 1
,∀ 1 ∈
0

Percebe-se que no segundo modo 2, que para qualquer valor de 1

com 2 = 0, a plataforma vibra com qualquer amplitude e o disco permanece

parado.

Para a análise gráfica dos resultados, foram dadas condições iniciais no

primeiro modo de vibração do sistema.
Tabela 2 – Condições Iniciais

Condições Iniciais
Deslocamento Inicial da Plataforma
1m
( = )
Velocidade Inicial da Plataforma
0 m/s
( = )
Deslocamento Inicial do Disco
2m
( = )
Velocidade Inicial do Disco ( = ) 0 m/s

1º Caso: Vibrações Livres sem Amortecimento e sem Forçamento

Tabela 3 – Parâmetros e Evolução – Vibrações Livres

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 0 N.s/m

Frequência de forçamento  0 rad/s

Amplitude de forçamento � 0N
Evolução
( ) 0 m (número de passos) 5000 ∆ ( ) 0,02
 Gráfico 1 – Deslocamento x Tempo – Vibrações Livres

Gráfico 2 - Deslocamento e Velocidade x Tempo (Plataforma) – Vibrações Livres

 Gráfico 3 - Deslocamento e Velocidade x Tempo (Disco) – Vibrações Livres

Como se observa no gráfico 1, dada a condição inicial no primeiro modo

o sistema só vibra neste primeiro modo, com a plataforma tendo amplitude
variando de -1m a 1m e o disco variando de -2m a 2m, com o sistema
oscilando em fase.
De acordo com os gráficos 1, 2 e 3, como era de se esperar, o sistema
vibra livremente e infinitamente, sem perdas de energia nem forçamentos. A
plataforma e o disco oscilam em suas frequências naturais. Percebe-se que a
amplitude de velocidades do disco é maior do que o da plataforma neste modo
de vibração.
O gráfico 4 mostra a combinação dos gráficos 1, 2 e 3 simultaneamente.
 Gráfico 4 – Deslocamento x Velocidade x Tempo – Vibrações Livres

O gráfico 5 abaixo mostra pontos importantes do sistema. Para os dois

corpos percebe-se que para o deslocamento zero (posição de equilíbrio) a
velocidade dos dois corpos é máxima e para o deslocamento máximo, 1m para
a plataforma e 2m para o disco, a velocidade é nula, ponto no qual as
acelerações são máximas. Este comportamento representa o comportamento
de um sistema oscilatório livre.

Gráfico 5 – Velocidade x Deslocamento – Vibrações Livres

2º Caso: Vibrações Livres Amortecidas sem Forçamento
Tabela 4 – Parâmetros e Evolução – Vibrações Livres Amortecidas

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 1 N.s/m

Frequência de forçamento  0 rad/s

Amplitude de forçamento � 0N
Evolução
( ) 0 m (número de passos) 5000 ∆ ( ) 0,02

Gráfico 6 – Deslocamento x Tempo – Vibrações Livres Amortecidas

Gráfico 7 – Deslocamento e Velocidade x Tempo (Plataforma) – Vibrações Livres
Amortecidas

Gráfico 8 – Deslocamento e Velocidade x Tempo (Disco) – Vibrações Livres

Amortecidas
 De acordo com os gráficos 6, 7 e 8, o comportamento agora com
amortecimento é bem diferente daquele de vibrações sem amortecimento.
Percebe-se que a adição de amortecimento no sistema provoca um
decaimento exponencial tanto no deslocamento quanto na velocidade,
tendendo a zero (até parar).
Nesse sistema, agora amortecido e sem forçamento, tanto a plataforma
quanto o disco oscilam, mas perdem energia. Uma vez iniciado o movimento,
os corpos não atingem mais seus deslocamentos máximos, onde os mesmos
tendem a parar.
Verifica-se também que para = 1 o sistema é sub-amortecido com
comportamento periódico (0 <  < 1).

O gráfico 9 mostra a combinação dos gráficos 1, 2 e 3 simultaneamente.

Gráfico 9 - Deslocamento x Velocidade x Tempo – Vibrações Livres

 Gráfico 10 -– Velocidade x Deslocamento – Vibrações Livres Amortecidas

Observa-se pelo gráfico 10 que com a presença de amortecimento há

perda de energia, o gráfico tem a forma de espiral, onde observa-se a perda de
velocidade tanto da plataforma quanto do disco.
Pelo gráfico, os corpos tendem a posição de equilíbrio, onde na forma de
espiral, as velocidades tendem a zero a partir da amplitude máxima, 1m para a
plataforma e 2m para o disco.

3º Caso: Vibrações Amortecidas com Forçamento Harmônico

Tabela 4 – Parâmetros e Evolução – Vibrações Livres Amortecidas

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 1 N.s/m

Frequência de forçamento  1 rad/s

Amplitude de forçamento � 1N
 Evolução
( ) 0 m (número de passos) 5000 ∆ ( ) 0,03

Gráfico 11 - Deslocamento x Tempo – Vibrações Amortecidas com Forçamento

Harmônico
 Gráfico 12 – Velocidade x Tempo (Plataforma) – Vibrações Amortecidas com
Forçamento Harmônico

Gráfico 13 – Velocidade x Tempo (Disco) – Vibrações Amortecidas com Forçamento

Harmônico

Pelos gráficos 11, 12 e 13, se observa o comportamento do

deslocamento e velocidade da plataforma e do disco diante forçamento
harmônico e amortecimento.
O sistema, oscila e perde energia devido ao amortecimento, mas como
há a presença de um forçamento externo, ao contrário das vibrações livres
amortecidas, o sistema permanece oscilando com amplitude constante, ou
seja, há um decaimento da amplitude dos deslocamentos e velocidades
durante um período de tempo e depois permanecem constantes.
Os deslocamentos ( ) tendem a solução particular ( ), onde a
amplitude passa a ser constante em regime permanente. O movimento do
sistema começa em regime transiente, mas em um determinado período de
tempo passa a ser permanente devido ao forçamento.

4º Caso: Ressonância
 O fenômeno da ressonância ocorre quando se tem um forçamento
harmônico, cuja freqüência de forçamento () coincide com uma das
frequências naturais do sistema ( ). Em razão disso o sinal ( ) de
deslocamento é amplificado, no caso sem amortecimento ele cresce
indefinidamente e com amortecimento o crescimento é limitado.
O estudo deste fenômeno é muito importante, pois muitos sistemas
vibratórios, não só mecânicos, entram em colapso devido ao fenômeno de
ressonância.
Como calculada anteriormente, a frequência fundamental do sistema é:

1 =
5
� / = 0,447213595 � /
5

Logo, a mesma será utilizada para o estudo da ressonância do sistema.

Assim, para colocar o sistema em ressonância deve-se igualar a freqüência de
forçamento () com a frequência natural (1 ).

4.1 Ressonância em Sistema Não Amortecido

Tabela 5 – Parâmetros e Evolução – Ressonância Não Amortecido

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 0 N.s/m

Frequência de forçamento  0,447213595 rad/s

Amplitude de forçamento � 1N
Evolução
( ) 0 m (número de passos) 10000 ∆ ( ) 0,03
 Gráfico 14 – Deslocamento x Tempo – Ressonância Não Amortecido

Gráfico 15 – Velocidade x Tempo (Plataforma) – Ressonância Não Amortecido

 Gráfico 16 – Velocidade x Tempo (Disco) – Ressonância Não Amortecido

Analisando os gráficos 14, 15 e 16, percebe-se que ao igualar a

freqüência de forçamento com a frequência fundamental, o sistema entra em
ressonância.
Percebe-se, que por se tratar de um sistema não amortecido, a
amplificação dos deslocamentos e velocidades da plataforma e do disco cresce
indefinidamente. As respostas ( ) em sua amplitude crescem de forma
ilimitada.

O gráfico 17 combina os gráficos 14, 15 e 16:

 Gráfico 17 - Deslocamento x Velocidade x Tempo – Ressonância Não
Amortecido

4.2 Ressonância em Sistema Amortecido

Tabela 6 – Parâmetros e Evolução – Ressonância Amortecido

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 1 N.s/m

Frequência de forçamento  0,447213595 rad/s

Amplitude de forçamento � 1N
Evolução
( ) 0 m (número de passos) 10000 ∆ ( ) 0,03
 Gráfico 18 – Deslocamento x Tempo – Ressonância Amortecido

Gráfico 19 – Velocidade x Tempo (Plataforma) – Ressonância Amortecido

 Gráfico 20 – Velocidade x Tempo (Disco) – Ressonância Amortecido

Pela análise dos gráficos 17, 18 e 19, diferentemente da ressonância em

sistema sem amortecimento, neste caso, com amortecimento, a resposta ( )
aumenta em sua amplitude de forma limitada, mantendo-se constante. A
amplificação das velocidades e deslocamentos do disco e plataforma cresce de
forma limitada. O gráfico 21 combina os gráficos 17, 18 e 19:
 Gráfico 21 – Velocidade x Deslocamento x Tempo – Ressonância Amortecido

Para a freqüência 2 = 1 � / verifica-se um fato curioso. Ao igualá-la

à frequência de forçamento, somente há crescimento da amplitude da
plataforma, mantendo-se o disco com amplitude constante. Analisando o
fenômento para o caso sem amortecimento, tem-se:

Tabela 7 - Parâmetros e Evolução – Ressonância em 2

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 0 N.s/m

Frequência de forçamento  1 rad/s

Amplitude de forçamento � 1N
Evolução
( ) 0 m (número de passos) 10000 ∆ ( ) 0,03

Gráfico 22 – Velocidade x Tempo – Ressonância em 2

Gráfico 23 – Velocidade x Tempo (Plataforma) – Ressonância em 2

Gráfico 24 - Velocidade x Tempo (Disco) – Ressonância em 2

 Analisando os gráficos 22, 23 e 24, percebe-se que nesta freqüência a
ressonância ocorre de forma a amplificar indefinidamente o deslocamento e
velocidades da plataforma. O disco por sua vez não tem amplificação
significativa, mantendo-se praticamente constante em termos de velocidade e
deslocamento.

5º Caso – Batimento

O fenômeno do batimento ocorre quando a freqüência de forçamento

externo () é muito próxima de uma das freqüências naturais do sistema ( ).
Considerando a frequência fundamental do sistema para a análise do
batimento,

1 =
5
� / = 0,447213595 � /
5
Adotando uma freqüência de forçamento bem próxima de 1 , analisa-se
o batimento para o caso sem amortecimento (regime permanente).
Tabela 8 - Parâmetros e Evolução – Batimento

Parâmetros do Sistema Mecânico

Massa da plataforma � 1 kg

Massa do disco � 2 kg

Rigidez das molas � 1 N/m

Constante de amortecimento � 0 N.s/m

Frequência de forçamento  0,42 rad/s

Amplitude de forçamento � 1N
Evolução
( ) 0 m (número de passos) 10000 ∆ ( ) 0,05
 Gráfico 25 – Deslocamento x Tempo (Plataforma) – Batimento

Gráfico 26 – Deslocamento x Tempo (Disco) - Batimento

Os gráficos 25 e 26 representam o fenômeno batimento, que se

assemelha ao “batimento de um coração”.
 Conclusão

Conclui-se que a aplicação do método numérico RK-4, diante dos

resultados encontrados, está de acordo com o esperado com a teoria. De
acordo com as análises feitas utilizando os princípios físicos e teóricos, os
resultados foram satisfatórios, uma vez que mostraram fielmente os estudos
feitos em sala de aula.
O método Runge-Kutta 4ª Ordem se mostrou muito preciso nas análises.
O estudo do sistema massa-mola-amortecedor mostrou varias situações que
poderiam ocorrer no mesmo, como presença ou não de amortecimento e
forçamento, ressonância e batimento.
Foram determinadas teoricamente as frequências naturais e modos de
vibração utilizando autovalores e autovetores e depois os mesmos foram
aplicados numericamente para ver se eram satisfeitos, e de acordo com os
resultados numéricos, estão de acordo com a teoria. A aplicação de C.I’s no
modo de vibração fundamental verificou que o calculado pelo RK-4 está de
acordo com o esperado, já que as condições iniciais nos modos fizeram o
sistema só oscilar nos mesmos.
Com as frequencias calculadas teoricamente, verificou-se
numericamente o fenômeno de ressonância e batimento do sistema, que
estavam de acordo com a fundamentação teórica.

Referências Bibliográficas

http://monografias.poli.ufrj.br/monografias/monopoli10005569.pdf

INMAN, D. J., Engeneering Vibration, 3th ed, New Jersey: Pearson

Education,Inc.,Upper Saddle River, 2008.
Anexos

“Script” Principal RK4SYSMain.m

% RK4 p/ sistemas de EDO

clear all; clc; close all

a = 0; % tempo inicial (s)

b = 100; % tempo final (s)

m = 5000; % número de intervalos

ya(1) = 1; % condição inicial de deslocamento da plataforma [m]

ya(2) = 0; % condição inicial de velocidade da plataforma [m/s]

ya(3) = 2; % condição inicial de deslocamento do disco [m]

ya(4) = 0; % condição inicial de velocidade do disco [m/s]

[X,Y] = RK4SYS ('secondode', a, b, ya, m);

% Plotagem do Resultado 3D

figure('Name','Gráfico 5');

plot3(X,Y(1,:),Y(2,:),'r', X, Y(3,:),Y(4,:),'b');grid;

xlabel ('Tempo [s]','fontsize',14);

ylabel ('Deslocamento [m]','fontsize',14);

zlabel ('Velocidade [m/s]','fontsize',14);

title('Deslocamento X Velocidade X Tempo ','fontsize',16);

legend('Plataforma','Disco');

% Plotagem do Resultado Deslocamento e Velocidade x Tempo do Disco

figure('Name','Gráfico 4');

plot(X,Y(3,:), 'r',X,Y(4,:), 'g');

xlabel ('Tempo(s)','fontsize',14);

ylabel ('Deslocamento [m], Velocidade [m/s]','fontsize',14);

title('Deslocamento e Velocidade X Tempo do Disco','fontsize',16);

legend('Deslocamento [m]','Velocidade [m/s]');

% Plotagem do Resultado Deslocamento e Velocidade X Tempo da Plataforma

figure('Name','Gráfico 3');

plot(X,Y(1,:),'r',X,Y(2,:), 'b');

xlabel ('Tempo(s)','fontsize',14);

ylabel ('Deslocamento [m], Velocidade [m/s]','fontsize',14);

title('Deslocamento e Velocidade X Tempo da Plataforma','fontsize',16);

legend('Deslocamento [m]','Velocidade [m/s]');

% Plotagem do Resultado Velocidade x Deslocamento (Plataforma e Disco)

figure('Name','Gráfico 2');

plot(Y(1,:),Y(2,:), 'r',Y(3,:),Y(4,:),'b');
xlabel ('Deslocamento [m]','fontsize',14);

ylabel ('Velocidade [m/s]','fontsize',14);

title('Velocidade X Deslocamento','fontsize',16);

legend('Plataforma','Disco');

% Plotagem do Resultado Deslocamento x Tempo (Plataforma e Disco)

figure('Name','Gráfico 1');

plot(X,Y(1,:),'r',X,Y(3,:),'b');

xlabel ('Tempo [s]','fontsize',14);

ylabel ('Deslocamento [m]','fontsize',14);

title('Deslocamento x Tempo ','fontsize',16);

legend('Plataforma','Disco');

“Script” Função do Runge-Kutta 4 RK4SYS.m

% Função do Runge-Kutta

function [X,Y] = RK4SYS(secondode, a, b, ya, m)

h = (b-a)/m; % calculo do passo (s)

N = size(ya, 2); % N recebe o numero de ya's na segunda ordem

X = zeros(1, m+1); % zera a matriz de tempo

Y = zeros(N, m+1); % zera a matriz de variaveis de estado

for K = 1:N

Y(K,1) = ya(K);

end

for J = 1:m

xj = X(J);

yj(1:N) = Y(1:N, J);

K1(1:N) = h*feval(secondode, xj, yj);

K2(1:N) = h*feval(secondode, xj+h/2, yj+K1/2);

K3(1:N) = h*feval(secondode, xj+h/2, yj+K2/2);

K4(1:N) = h*feval(secondode, xj+h, yj+K3);

Y(1:N, J+1) = yj(1:N) + (K1(1:N) + 2*K2(1:N) + 2*K3(1:N) + K4(1:N))/6;

X(J+1) = a+h*J;

end

return

“Script” das EDO’s secondode.m

function [fz] = secondode(x, y)

M = 1; % massa do corpo 1 (plataforma) [Kg]

m = 2; % massa do corpo 2 (disco) [Kg]

k = 1; % rigidez da mola [N/m]

c = 1; % constante de amortecimento [N.s/m]

Omega = 1; % frequencia de forçamento [rad/s]

F = 1; % amplitude de forçamento [N]

% sistema de equações

fz(1) = y(2);

fz(2) = ((k*(2*y(3)-5*y(1)))/(3*M+m))+((2*c*(y(4)-
y(2)))/(3*M+m))+(3*F*cos(Omega*x)/(3*M+m));

fz(3) = y(4);

fz(4) = (k*(2*M-m)*y(1)-2*M*y(3)+c*2*M*(y(2)-
y(4))+(m*F*cos(Omega*x)))/(m*(3*M+m));

return;