#!/usr/bin/env python
# coding: utf-8

# Universidade Federal do Rio Grande do Sul (UFRGS)   
# Programa de Pós-Graduação em Engenharia Civil (PPGEC)   
# 
# # PEC00025: Introduction to Vibration Theory
# 
# 
# ### Class 18 - Lagrange's equation
# 
# [1.   The Lagrange's equation](#section_1)  
# [2.   Example: generic s.d.o.f. system](#section_2)  
# [3.   Example: generic m.d.o.f. system](#section_3)  
# [4.   Example: simply supported beam](#section_4)  
# [4.1. Discrete flexibility and distributed mass](#section_41)  
# [4.2. Discrete mass and distributed flexibility](#section_42)  
# [4.3. Distributed mass and distributed flexibility](#section_43)  
# [5.   Example: nonlinear simple pendulum](#section_5)  
# [6.   Example: tower with flexible foundation](#section_6)  
# [7.   Example: nonlinear double pendulum](#section_7)  
# [8.   Assignments](#section_8)  
# 
# ---
# _Prof. Marcelo M. Rocha, Dr.techn._ [(ORCID)](https://orcid.org/0000-0001-5640-1020)  
# _Porto Alegre, RS, Brazil_ 
# 

# In[41]:


# Importing Python modules required for this notebook
# (this cell must be executed with "shift+enter" before any other Python cell)

import numpy as np
import matplotlib.pyplot as plt
from   scipy.integrate import odeint


# ## 1. The Lagrange's equation <a name="section_1"></a> 
# 

# In[ ]:





# ## 2. Example: generic s.d.o.f. system <a name="section_2"></a> 
# 

# In[ ]:





# ## 3. Example: generic m.d.o.f. system <a name="section_3"></a> 
# 

# In[ ]:





# ## 4. Example: simply supported beam <a name="section_4"></a> 
# 
# Now we shall use the Lagrange's equation for solving a simply supported beam, considering
# three different models:
# 
# 1. Distributed mass and discrete flexibility,
# 2. Distributed flexibility and discrete mass, and
# 3. Distributed mass and distributed flexibility.
# 
# <img src="images/simply_supported_lagrange.png" alt="Simply supported beam" width="360px"/>
# 
# The three beam models are useful for illustrating the use of Lagrange's equation with
# different energy evaluation approaches. At the end, the three solutions may be
# compared for evaluating the best estimates of discrete parameters $k_\theta$ and $M$
# such that the three models give the same beam natural vibration frequency.
# 
# 
# ### 4.1. Distributed mass and discrete flexibility <a name="section_41"></a> 
# 
# For the sake of compatibility with the other two models, we define the generalized coordinate
# $u$ as being the vertical displacement at beam center. The potential elastic energy
# depends on the relative rotation, $\theta$, between the two segments:
# 
# $$ \theta = \arcsin \left( \frac{2u}{L}\right) \approx \frac{2u}{L}$$
# 
# where the small displacements assumption ($\sin \beta \approx \beta$) has been used. Hence:
# 
# $$ V = \frac{1}{2} k_\theta \theta^2 \approx 
#                  \frac{2 k_\theta}{L^2} u^2$$
# 
# The kinetic energy is obtained by integrating the velocities along the two bars:
# 
# $$ T = 2 \cdot \frac{1}{2} \int_0^{L/2} {\left( \dot{r}_x^2 + \dot{r}_y^2 \right) \; \mu d\xi} $$
# 
# with:
# 
# \begin{align*}
#  r_x &= \xi \cos\left( \frac{2u}{L}\right) \\
#  r_y &= \xi \sin\left( \frac{2u}{L}\right)
# \end{align*}
# 
# Observe that, without the small displacements assumption, the dummy variable $\xi$ 
# is a line differential, not the projection in $x$ direction.
# Corresponding time derivatives are:
# 
# \begin{align*}
# \dot{r}_x &= \frac{dr_x}{du} \; \frac{du}{dt} =
#             -\frac{2\xi}{L} \sin\left( \frac{2u}{L}\right) \, \dot{u}
#     \approx -\frac{4 u \dot{u}}{L^2} \xi \approx 0 \\
# \dot{r}_y &= \frac{dr_y}{du} \; \frac{du}{dt} = 
#              \frac{2\xi}{L} \cos\left( \frac{2u}{L}\right) \, \dot{u}
#     \approx  \frac{2\dot{u}}{L} \xi \\
# \end{align*}
# 
# where the small displacements assumption is made ($\sin \beta \approx \beta$ and
# $\cos \beta \approx 1$) and the velocity in direction $x$ is
# neglected. Replacing these results in the kinetic energy yields:
# 
# $$ T \approx \frac{4\dot{u}^2}{L^2} \int_0^{L/2} {\xi^2 \; \mu d\xi}
#            = \frac{\mu L}{6} \, \dot{u}^2 $$
# 
# Now the two terms of Lagrange's equation can be calculated. Firstly the kinectic energy:
# 
# $$ \frac{\partial T}{\partial \dot{u}} = \frac{\mu L}{3} \, \dot{u} $$
# 
# and then for the potential elastic energy:
# 
# $$ \frac{\partial V}{\partial u} = \frac{4 k_\theta }{L^2} \, u $$
# 
# Applying Lagrange's equation gives finally the equilibrium equation:
# 
# $$ \frac{\mu L}{3} \, \ddot{u} +  \frac{4 k_\theta }{L^2} \, u = 0 $$
# 
# which means that the natural vibration frequency is:
# 
# $$ \omega_{\rm n} = \sqrt{\frac{12 k_\theta}{\mu L^3}} $$
# 

# 
# ### 4.2. Distributed flexibility and discrete mass <a name="section_42"></a> 
# 
# In this case we must assume a deflected shape. A good option is the elastic line resulting
# from a point load at beam center:
# 
# $$ w(\xi) = \frac{P L^3}{12 EI} \left( \frac{3}{4} \xi - \xi^3 \right), 
#             \hspace{1cm} \xi \leq \frac{1}{2} $$
#             
# where the variable $\xi$ is a non-dimensional coordinate, $\xi = x/L$.
# This function is valid up to beam center and its symmetry must be accounted for. 
# The assumed shape must be rescaled, such that the displacement at beam center is our generalized
# coordinate, $u$. The maximumum central deflection from the function above is:
# 
# $$ w_{\rm max} = \frac{P L^3}{48 EI} $$
# 
# Re-scaling gives:
# 
# $$ w(\xi) = 4 u \left( \frac{3}{4} \xi - \xi^3 \right), 
#             \hspace{1cm} \xi \leq \frac{1}{2} $$
# 
# The rotation and the curvature functions are, respectively:
# 
# \begin{align*}
#  w^{\prime}(\xi)       &=  \frac{ 4}{L  } u \left( \frac{3}{4} - 3\xi^2 \right) \\
#  w^{\prime\prime}(\xi) &= -\frac{24}{L^2} u \xi 
# \end{align*}
# 
# where one must be aware that:
# 
# $$ \frac{d}{dx} = \frac{d}{d\xi}\,\frac{d\xi}{dx} = \frac{1}{L} \, \frac{d}{d\xi}$$
# 
# The potential elastic energy is then given by:
# 
# $$ V = 2 \cdot \frac{1}{2} \int_{0}^{1/2} {EI \left[ w^{\prime\prime}(\xi) \right]^2 \, L d\xi}
#      = \frac{24 EI}{L^3} \, u^2 $$
# 
# On the other hand, the kinetic energy evaluation is straighforward:
# 
# $$ T = \frac{1}{2} M\dot{u}^2 $$
# 
# Now we calculate the two terms of Lagrange's equation, starting with the kinectic energy:
# 
# $$ \frac{\partial T}{\partial \dot{u}} = M \, \dot{u} $$
# 
# and then for the potential elastic energy:
# 
# $$ \frac{\partial V}{\partial u} = \frac{48 EI}{L^3} \, u $$
# 
# Applying Lagrange's equation gives finally the equilibrium equation:
# 
# $$ M \, \ddot{u} +  \frac{48 EI}{L^3} \, u = 0 $$
# 
# which means that the natural vibration frequency is:
# 
# $$ \omega_{\rm n} = \sqrt{\frac{48 EI}{ML^3}} $$
# 

# 
# ### 4.3. Distributed mass and distributed flexibility <a name="section_43"></a> 
# 
# Here we assume a half-sine as the deflected shape (that we _know_ to be the right solution):
# 
# $$ w(\xi) = u \sin\left( \pi \xi \right) $$
#             
# The rotation and the curvature functions are, respectively:
# 
# \begin{align*}
#  w^{\prime}(\xi)       &=  \frac{\pi   u}{L  } \, \cos \left(\pi \xi \right) \\
#  w^{\prime\prime}(\xi) &= -\frac{\pi^2 u}{L^2} \, \sin \left(\pi \xi \right)
# \end{align*}
# 
# The velocity is assumed (small displacements) to have only the vertical component,
# which is then calculated as:
# 
# $$ \dot{w}(\xi) = \frac{dw}{du} \, \frac{du}{dt} = \dot{u} \sin\left( \pi \xi \right)$$
# 
# The potential elastic energy is then given by:
# 
# $$ V = \frac{1}{2} \int_{0}^{1} {EI \left[ w^{\prime\prime}(\xi) \right]^2 \; L d\xi}
#      = \frac{\pi^4 EI}{4 L^3} \, u^2 $$
# 
# A comparison with the previous model leads to $\pi^4/4 \approx 24$, such
# that the potential energy for both models are quite close (error is $\approx 1.5$%)
# even by using different functions for the deflected shape. 
# Comparison with the discrete flexibility model demands that:
# 
# $$ \frac{2 k_\theta}{L^2} = \frac{\pi^4 EI}{4 L^3} $$
# 
# and consequently the hinge flexibility to match the continuous beam flexibility
# must be:
# 
# $$ k_\theta = \frac{\pi^4 EI}{8 L}$$
# 
# On the other hand, the kinetic energy is given by:
# 
# $$ T = \frac{1}{2} \int_{0}^{1} {\mu \left[ \dot{w}(\xi) \right]^2 \; L d\xi}
#      = \frac{\mu L}{4} \, \dot{u}^2 $$
# 
# A comparison with the discrete mass model demands that:
# 
# $$ \frac{M}{2} = \frac{\mu L}{4}$$
# 
# and consequently the discrete mass must be $M = \mu L/2$, which is half the total beam mass.
# 
# Finally, taking derivatives and applying Lagrange's equation gives the equilibrium equation:
# 
# $$ \frac{\mu L}{2} \, \ddot{u} +  \frac{\pi^4 EI}{2 L^3} \, u = 0 $$
# 
# and the natural vibration frequency results as :
# 
# $$ \omega_{\rm n} = \left( \frac{\pi}{L} \right)^2 \sqrt{\frac{EI}{\mu}} $$
# 
# which is the _exact_ solution for the first natural frequency of a simply supported beam.
# 

# ## 5. Example: simple pendulum <a name="section_5"></a> 
# 
# 
# <img src="images/simple_pendulum.png" alt="Simple pendulum" width="240px"/>
# 
# 

# In[ ]:





# ## 6. Example: tower with flexible foundation <a name="section_6"></a> 
# 

# In[ ]:





# ## 7. Example: nonlinear double pendulum <a name="section_7"></a> 
# 
# 
# <img src="images/double_pendulum.png" alt="Double pendulum" width="240px"/>
# 
# 
# Define $d\vec{y}/dt$ function:
# 

# In[42]:


def double_pendulum(y, t, m1, h1, m2, h2, g):
    
    th1, w1, th2, w2 = y
    
    mt  = m1  + m2
    dth = th1 - th2
    
    cd  = np.cos(dth)
    sd  = np.sin(dth)
    s1  = np.sin(th1)
    s2  = np.sin(th2)
    
    A   = np.array([[mt*h1,    m2*h2*cd], 
                    [m2*h1*cd, m2*h2  ]])

    B   = np.array([-m2*h2*w2*w2*sd - g*mt*s1,
                     m2*h1*w1*w1*sd - g*m2*s2])
    
    a   = np.linalg.solve(A, B)
    
    return [w1, a[0], w2, a[1]]


# Call ```scipy``` for differential equation solution:
# 

# In[49]:


t  =  np.linspace(0, 4, 16384)

m1 =  1.0;   h1 = 1.0
m2 =  1.5;   h2 = 0.5

g  =  9.81

Y0 = [np.pi+0.001, 0., np.pi/2, 0.]

plt.figure(1, figsize=(12, 5.5), clear=True)

for k in range(2):
    
    Y  =  odeint(double_pendulum, Y0, t, args=(m1, h1, m2, h2, g))

    x1 =  h1*np.sin(Y[:,0])
    y1 = -h1*np.cos(Y[:,0])

    x2 =  x1 + h2*np.sin(Y[:,1])
    y2 =  y1 - h2*np.cos(Y[:,1])

    plt.subplot(1,2,k+1)
    plt.plot(x1,y1, 'r:')
    plt.plot(x2,y2, 'b', linewidth=0.5)

    plt.xlim(-2, 2);   plt.xlabel('X2 (m)') 
    plt.ylim(-2, 2);   plt.ylabel('Y2 (m)') 
    plt.title('Simulation {0}'.format(k+1))

    plt.grid(True) 
    
    Y0 = [np.pi+0.002, 0, np.pi/2, 0.]   # repeat with butterfly wing flap


# ## 8. Assignments <a name="section_8"></a> 
# 

# In[ ]: