CIE_4140_Lecture_2_1_Python

CIE_4140_Lecture_2_1_Python#

import sympy as sp
u = sp.symbols('u',cls=sp.Function)
t, omega_0 = sp.symbols('t omega_0')
C, s = sp.symbols('C, s')

Equation governing free motion of an SDOF

Equation_of_Motion= sp.diff(u(t),t,2)+omega_0 **2 * u(t)
display(Equation_of_Motion)
\[\displaystyle \omega_{0}^{2} u{\left(t \right)} + \frac{d^{2}}{d t^{2}} u{\left(t \right)}\]

Characteristic equation

CharacteristicEquation = sp.simplify(Equation_of_Motion.subs(u(t), C*sp.exp(s*t)) / (C * sp.exp(s*t)))
display(CharacteristicEquation)
\[\displaystyle \omega_{0}^{2} + s^{2}\]

The Eigenvalues

EigenValues = sp.solve(CharacteristicEquation,s)
display(EigenValues)
[-I*omega_0, I*omega_0]

The Eigenfrequencies

EigenFrequency = []
for k in [0,1]:
    EigenFrequency.append(EigenValues[k]/ sp.I)
    display(EigenFrequency[k])
\[\displaystyle - \omega_{0}\]
\[\displaystyle \omega_{0}\]

The General Solution

C_1, C_2 = sp.symbols('C_1, C_2')
General_Solution_Complex = C_1*sp.exp(EigenValues[0]*t) + C_2*sp.exp(EigenValues[0]*t)
display(General_Solution_Complex)
\[\displaystyle C_{1} e^{- i \omega_{0} t} + C_{2} e^{- i \omega_{0} t}\]
A_1, A_2 = sp.symbols('A_1, A_2')
General_Solution_Real= A_1*sp.sin(EigenFrequency[0]*t)+A_2*sp.cos(EigenFrequency[1]*t)
display(General_Solution_Real)
\[\displaystyle - A_{1} \sin{\left(\omega_{0} t \right)} + A_{2} \cos{\left(\omega_{0} t \right)}\]
A, phi = sp.symbols('A phi')
General_Solution_Real_Compact =A*sp.cos(EigenFrequency[1]*t+phi)
display(General_Solution_Real_Compact)
\[\displaystyle A \cos{\left(\omega_{0} t + \phi \right)}\]

Visualisation

sp.plot(General_Solution_Real_Compact.subs([(A,1),(phi,0),(omega_0, 1.3)]),(t,0,10));
../_images/ddf47fa7681f03b55a15aeacba5526a8c7ee41326311240705251e939fbb2023.png
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
General_Solution_Real_Compact_func = sp.lambdify((t,A,phi,omega_0),General_Solution_Real_Compact)
fig, ax = plt.subplots()
tdata = np.linspace(0,10,100)
line, = ax.plot([], [])
ax.set_xlim(0, 10)
ax.set_ylim(-1, 1)

def update(frame):
    ydata = General_Solution_Real_Compact_func(t=tdata,A=1,phi=frame,omega_0=1.3)
    ax.set_title("$\phi$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0, 2*np.pi, 100),interval = 30)
plt.show()
fig, ax = plt.subplots()
tdata = np.linspace(0,10,100)
line, = ax.plot([], [])
ax.set_xlim(0, 10)
ax.set_ylim(-2, 2)

def update(frame):
    ydata = General_Solution_Real_Compact_func(t=tdata,A=frame,phi=0,omega_0=1.3)
    ax.set_title("$A$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(1, 2, 100),interval = 30)
plt.show()
fig, ax = plt.subplots()
tdata = np.linspace(0,10,100)
line, = ax.plot([], [])
ax.set_xlim(0, 10)
ax.set_ylim(-1, 1)

def update(frame):
    ydata = General_Solution_Real_Compact_func(t=tdata,A=1,phi=0,omega_0=frame)
    ax.set_title("$\omega_0$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0.3, np.pi, 100),interval = 30)
plt.show()