Considerando o problema do pêndulo simples mostrado na Figura.
Sendo a matriz de transformação é dada por:
\begin{equation} \textbf{T}_{\theta} = \begin{bmatrix} c\theta & s\theta & 0 \\ -s\theta & c\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{equation}O vetor posição descrito na base $\textit{S}_{1}$ é dada por:
\begin{equation} _{S_1}^{S_1}\textbf{p}^P = \begin{Bmatrix} 0 \\ L \\ 0 \end{Bmatrix} \end{equation}Enquanto a aceleração do pêndulo pode ser escrita conforme segue.
\begin{equation} _{S_1}^{I}\textbf{a}^P = \; _{S_1}^{I}\textbf{a}^O + \; _{S_1}\textbf{a}^{P/S_1} + \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \big( _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \: _{S_1}\textbf{p}^{P/S_1} \big) + 2 \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; _{S_1}\textbf{v}^{P/S_1} + _{S_1}^{I}\boldsymbol{\alpha}^{S_1} \times \; _{S_1}\textbf{p}^{P/S_1} \end{equation}A partir de uma análise do movimento do sistema observa-se que $\textbf{a}^{O} = \textbf{a}^{P/S_1} = 2 \; ^{I}\boldsymbol{\omega}^{S_1} \times \textbf{v}^{P/S_1} = 0$ e com isso, tem-se:
\begin{equation} _{S_1}^{I}\boldsymbol{\alpha}^{S_1} \times \; _{S_1}\textbf{p}^{P/S_1} = \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \ddot{\theta} \\ 0 & L & 0 \\ \end{bmatrix} = \begin{Bmatrix} 0 \\ -L\ddot{\theta} \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; \big( _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \, _{S_1}\textbf{p}^{P/S_1} \big) = \; _{S_1}^{I}\boldsymbol{\omega}^{S_1} \times \; \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \dot{\theta} \\ 0 & L & 0 \\ \end{bmatrix} = \begin{bmatrix} \textbf{e}_1^{(1)} & \textbf{e}_2^{(1)} & \textbf{e}_3^{(1)} \\ 0 & 0 & \dot{\theta} \\ -L\dot{\theta} & 0 & 0 \\ \end{bmatrix} = \begin{Bmatrix} 0 \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} \end{equation}Portanto, a aceleração é dada por:
\begin{equation} _{S_1}^{I}\textbf{a}^P = \begin{Bmatrix} -L\ddot{\theta} \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} \end{equation}Avaliando as forças no pêndulo, considera-se uma força resultante $F$ composta pelo peso, $P$, uma força dissipativa viscosa linear, $F_D$, uma excitação externa, $F_E$ e uma força na haste, $F_H$, conforme mostrado na Figura anterior.
Com isso, tem-se:
\begin{equation} _{S_1}\textbf{P} = \textbf{T}_{\theta}\; _{I}\textbf{P} = \begin{bmatrix} c\theta & s\theta & 0 \\ -s\theta & c\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{Bmatrix} 0 \\ mg \\ 0 \end{Bmatrix} = \begin{Bmatrix} mgs\theta \\ mgc\theta \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_D = \begin{Bmatrix} CL\dot{\theta} \\ 0 \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_E = \begin{Bmatrix} -F(t) \\ 0 \\ 0 \end{Bmatrix} \end{equation}\begin{equation} _{S_1}\textbf{F}_H = \begin{Bmatrix} 0 \\ -T \\ 0 \end{Bmatrix} \end{equation}Desta forma, aplicando a segunda lei de Newton, $m \;_{S_1}^{I}\textbf{a}^P = \,_{S_1}\textbf{F}$, tem-se que:
\begin{equation} m \; _{S_1}\textbf{a}^P = m \begin{Bmatrix} -L\ddot{\theta} \\ -L\dot{\theta}^2 \\ 0 \end{Bmatrix} = \begin{Bmatrix} mgs\theta + CL\dot{\theta} - F(t) \\ mgc\theta - T \\ 0 \end{Bmatrix} \end{equation}O que resulta nas seguintes equações:
\begin{equation} mL\ddot{\theta} + CL\dot{\theta} + mgs\theta = F(t) \\ T = mL\dot{\theta}^2 + mgc\theta \end{equation}Note que a segunda esquação estabelece o valor da tensão na haste enquanto a primeira equação define o movimento do pêndulo que, trata-se de uma equação diferencial não-linear. Considerando pequenos ângulos, $\theta ≪ 1$, pode-se assumir que $s\theta \approx \theta$. Essa simplificação usual transforma a equação de movimento em uma equação diferencial linear, que representa um oscilador linear:
\begin{equation} mL\ddot{\theta} + CL\dot{\theta} + mg\theta = F(t) \end{equation}A fim de obter a solução numérica da equação de movimento utilizando o método de Runge-Kutta de 4ª ordem. Fazemos a seguinte mudança de variáveis: $x_1 = \theta$, $x_2 = \dot{\theta}$, permitindo que a equação de movimento não-linear e linear sejam escritas na forma canônica, conforme expresso a seguir.
$\textbf{Não-linear}$
\begin{equation} \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \frac{F(t)}{mL} - \frac{C}{m}x_2 - \frac{g}{L}s(x_1) \end{cases} \end{equation}$\textbf{Linear}$
\begin{equation} \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \frac{F(t)}{mL} - \frac{C}{m}x_2 - \frac{g}{L}x_1 \end{cases} \end{equation}onde $\textbf{x}$ é o vetor de variáveis de estado dado por,
\begin{equation} \textbf{x} = \begin{Bmatrix} x_1 \\ x_2 \\ \end{Bmatrix} = \begin{Bmatrix} \theta \\ \dot{\theta} \\ \end{Bmatrix} \end{equation}Em seguida vamos definir uma função em Python para as equações de movimento não-linear e linear.
%matplotlib inline
from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
from RK4 import integrate
#============================================================================#
# Right-hand side of simple pendulum differential equations #
#============================================================================#
def lin_sim_pen(x, t):
ff = np.zeros_like(x)
ff[0] = x[1]
ff[1] = amp*sin(omega*t) - ksi*x[1] -(G/L)*x[0]
return ff
def nonlin_sim_pen(x, t):
ff = np.zeros_like(x)
ff[0] = x[1]
ff[1] = amp*sin(omega*t) - ksi*x[1] -(G/L)*sin(x[0])
return ff
Vamos criar uma função para plotar os resultados obtidos.
def plots(y, X1, X2, label, ls, ax1, ax2):
t = y[:,0] # time array
theta = y[:,1] # angle array
w = y[:,2] # angular velocity array
theta = np.degrees(theta) # converts from rad to degree
w = np.degrees(w) # converts from rad to degree
#========================================================================#
# Plot results #
#========================================================================#
ax1.plot(t, theta, ls, label=label, alpha=1.0)
ax1.set_xlabel(' Tempo (s) ')
ax1.set_ylabel('$\\theta [\degree]$')
ax1.set_title(' Pêndulo Simples')
ax1.set_xlim(t[0], t[-1])
ax1.set_ylim(-1.5*np.max(theta), 1.5*np.max(theta))
ax1.legend()
#
# State space
#
ax2.plot(theta, w, ls, label=label)
ax2.set_xlabel('$\\theta [\degree]$')
ax2.set_ylabel('$\\omega [\degree/s]$')
ax2.set_title(' Pêndulo Simples')
ax2.set_xlim(1.5*np.min(theta), 1.5*np.max(theta))
ax2.legend()
Agora, podemos definir os parâmtros do problema e utlizar as funções de integração numérica pelo método de Runge-Kutta, implementadas anteriormente, para obter a resposta do sistema. Por fim plotamos os resultados.
Primeiramente, assumimos como condições iniciais, uma condição de pequeno deslocamento, $\theta_0 = -10°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.0 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -10.0 # initial angle (degrees)
w1 = 0.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
Conforme pode ser observado, para pequenos ângulos, a resposta do sistema são semelhantes para a equação de movimento linear e não-linear.
Por outro lado, para ângulos maiores, digamos $\theta_{0} = -120°$ e mantendo $\dot{\theta}_0 = 0°/s$, o pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.0 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -120.0 # initial angle (degrees)
w1 = 0.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
Considerando neste momento as seguintes condições iniciais, $\theta_{0} = -120°$ e $\dot{\theta}_0 = 200 °/s$, onde uma velocidade angular não nula é assumida, tem-se o comparativo mostrado a seguir.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.0 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -120.0 # initial angle (degrees)
w1 = 200.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-200.0, 400);
ax2.set_xlim(-200.0, 1000);
Em seguida, considera-se uma velocidade angular iniciais ainda maior, $\theta_0=500°/𝑠$, mantendo o mesmo valor de deslocamento angular inicial. Observa-se que o pêndulo linear tem um comportamento peculiar, que $\textbf{não}$ tem nenhum sentido físico para um sistema sem nenhum tipo de forçamento externo. Isso ocorre porque é assumido a hipótese de pequenas rotações que não condiz com a situação analisada. Em seguida, note o comportamento esperado para o pêndulo quando considera-se a não-linearidade.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.0 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -120.0 # initial angle (degrees)
w1 = 500.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-300.0, 600);
ax2.set_xlim(-300.0, 1000);
Considerando um amortecimento viscoso, $\xi = 0.08$, assumimos como condições iniciais, uma condição de pequeno deslocamento, $\theta_0 = -10°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s. Note que a resposta do sistema são semelhantes e apresentam decaimento exponencial com o decorrer do tempo.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.08 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -10.0 # initial angle (degrees)
w1 = 0.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
Considerando um amortecimento viscoso, $\xi = 0.08$ e assumindo um ângulo maior, $\theta_0 = -120°$ e velocidade angular de $\dot{\theta}_0 = 0°$/s. O pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.08 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -120.0 # initial angle (degrees)
w1 = 0.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
Considerando um amortecimento viscoso, $\xi = 0.2$ e assumindo um ângulo maior, $\theta_0 = -120°$ e velocidade angular de $\dot{\theta}_0 = 500°$/s. O pêndulo simples é não-linear e é possível observar a influência na resposta do sistema, conforme pode ser visualizado nos resultados que seguem.
#========================================================================#
# Inputs #
#========================================================================#
G = 9.8 # acceleration due to gravity, in m/s^2
L = 1.0 # length of pendulum in m
ksi = 0.2 # viscous damping ratio
amp = 0.0 # excitation amplitude
omega = 0.0 # excitation frequency
t0 = 0. # initial time
dt = 0.05 # time step
N = 1000 # number of time steps
th1 = -120.0 # initial angle (degrees)
w1 = 500.0 # initial angular velocity (degrees per second)
init_conditions = np.radians([th1, w1]) # converts to radians
plt.close('all')
fig = plt.figure(1, figsize=(14,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#========================================================================#
# Nonlinear Simple Pendulum Solution using RK4 #
#========================================================================#
y= integrate(t0, dt, N, init_conditions, nonlin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'não-linear','-r', ax1, ax2)
#========================================================================#
# Linear Simple Pendulum Solution using RK4 #
#========================================================================#
y=integrate(t0, dt, N, init_conditions, lin_sim_pen) # matrix with results
t = y[:,0] # time
theta = y[:,1] # angle
w = y[:,2] # angular velocity
# converts to cartesian coordenates
X1 = -L*sin(theta)
X2 = L*cos(theta)
plots(y, X1, X2, 'linear', '-b', ax1, ax2)
ax1.set_ylim(-300.0, 1800);
ax2.set_xlim(-200.0, 1800);