Vibrations exercices

Exercise1.py

import numpy as np
import sympy as sp

sin, cos = sp.sin, sp.cos

k1, k2 = sp.symbols("k1 k2", real=True, positive=True)
L1, L2 = sp.symbols("L1 L2", real=True, positive=True)
a, p, m= sp.symbols("a p m", real=True, positive=True)
theta = sp.symbols("theta")
dtheta = sp.symbols("dtheta")


C0 = np.array([0, p, 0])
P0 = C0 + np.array([0, a, 0])
P1 = P0 - np.array([L1, 0, 0])
P2 = P0 + np.array([L2, 0, 0])

C = np.array([p*theta, p, 0])
P = C + a*np.array([sin(theta), cos(theta), 0])

# Energia potencial elastica
x1 = sp.sqrt(sum((P-P1)**2)) - sp.sqrt(sum((P0-P1)**2))
print("x1 = ", x1)
x1ser = sp.series(x1, theta, n=2).removeO()
print("x1 series = ", x1ser)

x2 = sp.sqrt(sum((P-P2)**2)) - sp.sqrt(sum((P0-P2)**2))
print("x2 = ", x2)
x2ser = sp.series(x2, theta, n=2).removeO()
print("x2 series = ", x2ser)

V = (k1*x1ser**2+k2*x2ser**2)/2
V = sp.simplify(V)
print("V = ", V)

# Energia cinetica
Ec1 = (1/2)*m*p**2  * dtheta**2
Ec2 = (1/2)*(1/2)*m*p**2  * dtheta**2
T = Ec1 + Ec2
print("T = ", T)

# Lagrangian
t = sp.symbols("t")
Theta = sp.Function("Theta")(t)
dTheta = sp.diff(Theta, t)
L = T - V
L = L.subs(theta, Theta)
L = L.subs(dtheta, dTheta)

equation = sp.diff(sp.diff(L, dTheta), t) - sp.diff(L, Theta)
print("Final Equation = ")
print(equation)

Exercise2.py

import sympy as sp
import numpy as np

sin, cos, pi = sp.sin, sp.cos, sp.pi

h = sp.symbols("h", positive = True, real=True)
if True:
	R, e = sp.symbols("R e", positive=True, real=True)
	Re = R + sp.Rational(1, 2)*e
	Ri = R - sp.Rational(1, 2)*e
else:
	Ri, Re = sp.symbols("Ri Re", positive=True, real=True)
	R = sp.Rational(1, 2) * (Ri+Re)
	e = Re-Ri

r, t, z = sp.symbols("r t z")
x = r * cos(t)
y = r * sin(t)
p = np.array([x, y, z])
pip = np.inner(p, p)
pop = np.tensordot(p, p, axes=0)

rho = 1
m = sp.integrate(rho*r, (r, Ri, Re))
m = sp.integrate(m, (t, pi/2, 3*pi/2))
m = sp.integrate(m, (z, -h/2, h/2))
m = sp.simplify(m)

V = sp.integrate(r, (r, Ri, Re))
V = sp.integrate(V, (t, pi/2, 3*pi/2))
V = sp.integrate(V, (z, -h/2, h/2))
V = sp.simplify(V)

CM = sp.integrate(rho*sp.Matrix(p)*r, (r, Ri, Re))
CM = sp.integrate(CM, (t, pi/2, 3*pi/2))
CM = sp.integrate(CM, (z, -h/2, h/2))
CM /= m
CM = sp.simplify(CM)

II = pip * sp.eye(3) - sp.Matrix(pop)
II = sp.integrate(rho*II*r, (r, Ri, Re))
II = sp.integrate(II, (t, pi/2, 3*pi/2))
II = sp.integrate(II, (z, -h/2, h/2))
II /= m
II = sp.simplify(II)

print("V = ", V)
print("CM = ", CM)
print("II = ", II)

a = -CM[0]
print("a = ", a)
g = sp.symbols("g")
m = sp.symbols("m")
theta = sp.symbols("theta")
dtheta = sp.symbols("dtheta")

# Energia potencial
V = -m*g*a*cos(theta)
print("V = ")
print(V)

# Energia cinetica
Ec1 = sp.Rational(1, 2)*m*Re**2 *dtheta**2
Ec2 = sp.Rational(1, 2)*m*R**2 *(1+sp.Rational(1, 4)*e**2/R**2)*dtheta**2
T = Ec1 + Ec2
T = sp.simplify(T)
print("T = ")
print(T)

# Lagrangiano
L = T - V
Theta = sp.Function("Theta")(t)
dTheta = sp.diff(Theta, t)
L = L.subs(theta, Theta)
L = L.subs(dtheta, dTheta)
equation = sp.diff(sp.diff(L, dTheta), t) - sp.diff(L, Theta)
equation = sp.simplify(equation)
print("Final Equation = ")
print(equation)