CIE_4140_Lecture_3_3_Python

Example 1: A harmonic load

import sympy as sp

u = sp.symbols('u',cls=sp.Function)
omega_n, F_0, m, Omega, zeta = sp.symbols('omega_n, F_0, m, Omega, zeta',real=True,positive=True)
t = sp.symbols('t',real=True)
u_0, v_0 = sp.symbols('u_0, v_0')

Equation of moation:

Equation_of_Motion= sp.Eq(sp.diff(u(t),t,2)+2*zeta*omega_n*sp.diff(u(t),t)+omega_n **2 * u(t),F_0*sp.cos(Omega*t)/m)
display(Equation_of_Motion)

\[\displaystyle \omega_{n}^{2} u{\left(t \right)} + 2 \omega_{n} \zeta \frac{d}{d t} u{\left(t \right)} + \frac{d^{2}}{d t^{2}} u{\left(t \right)} = \frac{F_{0} \cos{\left(\Omega t \right)}}{m}\]

Load:

sp.plot(1*sp.sin(1.1*t)/1,(t,0,200));

Calculation of Duhamel integral (for some reason this cell only succeeds on the second try,this bug is reported sympy/sympy#24518)

omega_1, t_tilda = sp.symbols('omega_1, t_tilda',real=True,positive=True)
u_forced = sp.simplify(F_0/m/omega_1 *sp.integrate(sp.sin(Omega*t_tilda)*sp.exp(-zeta*omega_n*(t-t_tilda))*sp.sin(omega_1*(t-t_tilda)),(t_tilda,0,t)))
display(u_forced)

\[\displaystyle \frac{F_{0} \left(\Omega^{3} \sin{\left(\omega_{1} t \right)} - \Omega^{2} \omega_{1} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} - \Omega \omega_{1}^{2} \sin{\left(\omega_{1} t \right)} - 2 \Omega \omega_{1} \omega_{n} \zeta e^{\omega_{n} t \zeta} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{1} \omega_{n} \zeta \cos{\left(\omega_{1} t \right)} + \Omega \omega_{n}^{2} \zeta^{2} \sin{\left(\omega_{1} t \right)} + \omega_{1}^{3} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} + \omega_{1} \omega_{n}^{2} \zeta^{2} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{1} \left(\Omega^{4} - 2 \Omega^{2} \omega_{1}^{2} + 2 \Omega^{2} \omega_{n}^{2} \zeta^{2} + \omega_{1}^{4} + 2 \omega_{1}^{2} \omega_{n}^{2} \zeta^{2} + \omega_{n}^{4} \zeta^{4}\right)}\]

u_forced = u_forced.subs(omega_1,omega_n * sp.sqrt(1-zeta**2))
display(u_forced)

\[\displaystyle \frac{F_{0} \left(\Omega^{3} \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - \Omega^{2} \omega_{n} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} + \Omega \omega_{n}^{2} \zeta^{2} \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - 2 \Omega \omega_{n}^{2} \zeta \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{n}^{2} \zeta \sqrt{1 - \zeta^{2}} \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - \Omega \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right) \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + \omega_{n}^{3} \zeta^{2} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} + \omega_{n}^{3} \left(1 - \zeta^{2}\right)^{\frac{3}{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{n} \sqrt{1 - \zeta^{2}} \left(\Omega^{4} + 2 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right) + \omega_{n}^{4} \zeta^{4} + 2 \omega_{n}^{4} \zeta^{2} \cdot \left(1 - \zeta^{2}\right) + \omega_{n}^{4} \left(1 - \zeta^{2}\right)^{2}\right)}\]

response = (sp.exp(-zeta*omega_n*t)*(u_0*sp.cos(omega_n*t)+
                                    (v_0/(omega_n * sp.sqrt(1-zeta**2))+
                                     u_0*zeta*omega_n/(omega_n * sp.sqrt(1-zeta**2)))*sp.sin(omega_n*t))+
            u_forced)
display(response)

\[\displaystyle \frac{F_{0} \left(\Omega^{3} \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - \Omega^{2} \omega_{n} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} + \Omega \omega_{n}^{2} \zeta^{2} \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - 2 \Omega \omega_{n}^{2} \zeta \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{n}^{2} \zeta \sqrt{1 - \zeta^{2}} \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} - \Omega \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right) \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + \omega_{n}^{3} \zeta^{2} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)} + \omega_{n}^{3} \left(1 - \zeta^{2}\right)^{\frac{3}{2}} e^{\omega_{n} t \zeta} \sin{\left(\Omega t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{n} \sqrt{1 - \zeta^{2}} \left(\Omega^{4} + 2 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right) + \omega_{n}^{4} \zeta^{4} + 2 \omega_{n}^{4} \zeta^{2} \cdot \left(1 - \zeta^{2}\right) + \omega_{n}^{4} \left(1 - \zeta^{2}\right)^{2}\right)} + \left(u_{0} \cos{\left(\omega_{n} t \right)} + \left(\frac{u_{0} \zeta}{\sqrt{1 - \zeta^{2}}} + \frac{v_{0}}{\omega_{n} \sqrt{1 - \zeta^{2}}}\right) \sin{\left(\omega_{n} t \right)}\right) e^{- \omega_{n} t \zeta}\]

sp.plot(response.subs([(u_0,0.1),(v_0,0.1),(F_0,1),(omega_n,1),(zeta,0.05),(m,1),(Omega,1.2)]), (t , 0 , 200));

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

response_func = sp.lambdify((u_0,v_0,F_0,omega_n,zeta,m,Omega,t),response)

fig, ax = plt.subplots()
tdata = np.linspace(0,200,500)
line, = ax.plot([], [])
ax.set_xlim(0, 200)
ax.set_ylim(-5, 5)

def update(frame):
    ydata = response_func(u_0=1,v_0=0,m=1,omega_n=1,F_0=1,zeta=0.05,Omega=frame,t=tdata)
    ax.set_title("Displacement versus time for $\Omega$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(1.1,1.6,100),interval = 100)
plt.show()

fig, ax = plt.subplots()
tdata = np.linspace(0,200,500)
line, = ax.plot([], [])
ax.set_xlim(0, 200)
ax.set_ylim(-9, 9)

def update(frame):
    ydata = response_func(u_0=1,v_0=0,m=1,omega_n=1,F_0=1,zeta=frame,Omega=1.1,t=tdata)
    ax.set_title("Displacement versus time for $\zeta$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

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

Example 2: a suddenly applied constant load

Equation of motion:

t = sp.symbols('t',real=True)
Equation_of_Motion= sp.Eq(sp.diff(u(t),t,2)+2*zeta*omega_n*sp.diff(u(t),t)+omega_n **2 * u(t),F_0*sp.Heaviside(t)/m)
display(Equation_of_Motion)

\[\displaystyle \omega_{n}^{2} u{\left(t \right)} + 2 \omega_{n} \zeta \frac{d}{d t} u{\left(t \right)} + \frac{d^{2}}{d t^{2}} u{\left(t \right)} = \frac{F_{0} \theta\left(t\right)}{m}\]

Load

sp.plot(sp.Heaviside(t-0),(t,-50,300));

Calculation of Duhamel integral

u_forced = sp.simplify(F_0/m/omega_1 *sp.integrate(sp.exp(-zeta*omega_n*(t-t_tilda))*sp.sin(omega_1*(t-t_tilda)),(t_tilda,0,t)))
display(u_forced)

\[\displaystyle \frac{F_{0} \left(\omega_{1} e^{\omega_{n} t \zeta} - \omega_{1} \cos{\left(\omega_{1} t \right)} - \omega_{n} \zeta \sin{\left(\omega_{1} t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{1} \left(\omega_{1}^{2} + \omega_{n}^{2} \zeta^{2}\right)}\]

The total solution

response = (sp.exp(-zeta*omega_n*t)*(u_0*sp.cos(omega_n*t)+
                                    (v_0/(omega_n * sp.sqrt(1-zeta**2))+
                                     u_0*zeta*omega_n/(omega_n * sp.sqrt(1-zeta**2)))*sp.sin(omega_n*t))+
            u_forced.subs(omega_1,omega_n*sp.sqrt(1-zeta**2)))
display(response)

\[\displaystyle \frac{F_{0} \left(- \omega_{n} \zeta \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + \omega_{n} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} - \omega_{n} \sqrt{1 - \zeta^{2}} \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{n} \sqrt{1 - \zeta^{2}} \left(\omega_{n}^{2} \zeta^{2} + \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right)\right)} + \left(u_{0} \cos{\left(\omega_{n} t \right)} + \left(\frac{u_{0} \zeta}{\sqrt{1 - \zeta^{2}}} + \frac{v_{0}}{\omega_{n} \sqrt{1 - \zeta^{2}}}\right) \sin{\left(\omega_{n} t \right)}\right) e^{- \omega_{n} t \zeta}\]

p0 = sp.plotting.plot(response.subs([(u_0,0),(v_0,0),(F_0,1),(m,1),(omega_n,1),(zeta,0   )]), (t , 0 , 50),label='$\zeta$=0'    ,legend=True,show=False)
p1 = sp.plotting.plot(response.subs([(u_0,0),(v_0,0),(F_0,1),(m,1),(omega_n,1),(zeta,0.05)]), (t , 0 , 50),label='$\zeta$=0.05' ,show=False)
p2 = sp.plotting.plot(response.subs([(u_0,0),(v_0,0),(F_0,1),(m,1),(omega_n,1),(zeta,0.1 )]), (t , 0 , 50),label='$\zeta$=0.1'  ,show=False)
p0.append(p1[0])
p0.append(p2[0])
p0.show()

response_func = sp.lambdify((u_0,v_0,F_0,omega_n,zeta,m,t),response)

fig, ax = plt.subplots()
tdata = np.linspace(0,200,500)
line, = ax.plot([], [])
ax.set_xlim(0, 200)
ax.set_ylim(0,2)

def update(frame):
    ydata = response_func(u_0=0,v_0=0,m=1,omega_n=1,F_0=1,zeta=frame,t=tdata)
    ax.set_title("Displacement versus time for $\zeta$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0.,0.2,100),interval = 100)
plt.show()

Example 3: a suddenly applied constant load of finite duration

Equation of motion

dt = sp.symbols('dt',real=True,positive=True)
Equation_of_Motion= sp.Eq(sp.diff(u(t),t,2)+2*zeta*omega_n*sp.diff(u(t),t)+omega_n **2 * u(t),F_0*(sp.Heaviside(t)-sp.Heaviside(t-dt))/m)
display(Equation_of_Motion)

\[\displaystyle \omega_{n}^{2} u{\left(t \right)} + 2 \omega_{n} \zeta \frac{d}{d t} u{\left(t \right)} + \frac{d^{2}}{d t^{2}} u{\left(t \right)} = \frac{F_{0} \left(\theta\left(t\right) - \theta\left(- dt + t\right)\right)}{m}\]

Load:

sp.plot(sp.Heaviside(t)-sp.Heaviside(t-4*sp.pi),(t,-50,30));

Calculation of the Duhamel integral

if \(t<dt\) then:

u_forced_1 = sp.simplify(F_0/m/omega_1 *sp.integrate(sp.exp(-zeta*omega_n*(t-t_tilda))*sp.sin(omega_1*(t-t_tilda)),(t_tilda,0,t)))
display(u_forced_1)

\[\displaystyle \frac{F_{0} \left(\omega_{1} e^{\omega_{n} t \zeta} - \omega_{1} \cos{\left(\omega_{1} t \right)} - \omega_{n} \zeta \sin{\left(\omega_{1} t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{1} \left(\omega_{1}^{2} + \omega_{n}^{2} \zeta^{2}\right)}\]

if \(t>dt\) then:

u_forced_2 = sp.simplify(F_0/m/omega_1 *sp.integrate(sp.exp(-zeta*omega_n*(t-t_tilda))*sp.sin(omega_1*(t-t_tilda)),(t_tilda,0,dt)))
display(u_forced_2)

\[\displaystyle \frac{F_{0} \left(\omega_{1} e^{dt \omega_{n} \zeta} \cos{\left(\omega_{1} \left(dt - t\right) \right)} - \omega_{1} \cos{\left(\omega_{1} t \right)} - \omega_{n} \zeta e^{dt \omega_{n} \zeta} \sin{\left(\omega_{1} \left(dt - t\right) \right)} - \omega_{n} \zeta \sin{\left(\omega_{1} t \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{1} \left(\omega_{1}^{2} + \omega_{n}^{2} \zeta^{2}\right)}\]

The total solution:

response = (sp.exp(-zeta*omega_n*t)*(u_0*sp.cos(omega_n*t)+
                                    (v_0/(omega_n * sp.sqrt(1-zeta**2))+
                                     u_0*zeta*omega_n/(omega_n * sp.sqrt(1-zeta**2)))*sp.sin(omega_n*t))+
            sp.Piecewise((u_forced_1.subs(omega_1,omega_n*sp.sqrt(1-zeta**2)),t<dt),
                         (u_forced_2.subs(omega_1,omega_n*sp.sqrt(1-zeta**2)),True)))
display(response)

\[\begin{split}\displaystyle \left(u_{0} \cos{\left(\omega_{n} t \right)} + \left(\frac{u_{0} \zeta}{\sqrt{1 - \zeta^{2}}} + \frac{v_{0}}{\omega_{n} \sqrt{1 - \zeta^{2}}}\right) \sin{\left(\omega_{n} t \right)}\right) e^{- \omega_{n} t \zeta} + \begin{cases} \frac{F_{0} \left(- \omega_{n} \zeta \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + \omega_{n} \sqrt{1 - \zeta^{2}} e^{\omega_{n} t \zeta} - \omega_{n} \sqrt{1 - \zeta^{2}} \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{n} \sqrt{1 - \zeta^{2}} \left(\omega_{n}^{2} \zeta^{2} + \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right)\right)} & \text{for}\: dt > t \\\frac{F_{0} \left(- \omega_{n} \zeta e^{dt \omega_{n} \zeta} \sin{\left(\omega_{n} \sqrt{1 - \zeta^{2}} \left(dt - t\right) \right)} - \omega_{n} \zeta \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + \omega_{n} \sqrt{1 - \zeta^{2}} e^{dt \omega_{n} \zeta} \cos{\left(\omega_{n} \sqrt{1 - \zeta^{2}} \left(dt - t\right) \right)} - \omega_{n} \sqrt{1 - \zeta^{2}} \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}\right) e^{- \omega_{n} t \zeta}}{m \omega_{n} \sqrt{1 - \zeta^{2}} \left(\omega_{n}^{2} \zeta^{2} + \omega_{n}^{2} \cdot \left(1 - \zeta^{2}\right)\right)} & \text{otherwise} \end{cases}\end{split}\]

response_func = sp.lambdify((u_0,v_0,F_0,omega_n,zeta,m,t,dt),response)

fig, ax = plt.subplots()
tdata = np.linspace(0,50,500)
line, = ax.plot([], [])
ax.set_xlim(0, 50)
ax.set_ylim(-2,2)

def update(frame):
    ydata = response_func(u_0=0,v_0=0,m=1,omega_n=1,F_0=1,zeta=0.05,dt=frame,t=tdata)
    ax.set_title("Displacement versus time for $dt$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(3*np.pi,5*np.pi,100),interval = 100)
plt.show()

fig, ax = plt.subplots()
tdata = np.linspace(0,50,500)
line, = ax.plot([], [])
ax.set_xlim(0, 50)
ax.set_ylim(-2,2)
i=0
for frame in np.linspace(3*np.pi,5*np.pi,100):
    def update(frame):
        ydata = response_func(u_0=0,v_0=0,m=1,omega_n=1,F_0=1,zeta=0.05,dt=frame,t=tdata)
        ax.set_title("Displacement versus time for $dt$ = "+str(np.round(frame,2)))
        line.set_data(tdata, ydata)
    update(frame)
    plt.savefig(f'img_{i}.png',transparent = False,facecolor = 'white')
    i+=1

import imageio
frames = []
for t in range(100):
    image = imageio.v2.imread(f'img_{t}.png')
    frames.append(image)

imageio.mimsave('./example.gif', # output gif
                frames,          # array of input frames
                fps = 10)         # optional: frames per second

fig, ax = plt.subplots()
tdata = np.linspace(0,50,500)
line, = ax.plot([], [])
ax.set_xlim(0, 50)
ax.set_ylim(-2,2)

def update(frame):
    ydata = response_func(u_0=0,v_0=0,m=1,omega_n=1,F_0=1,zeta=0,dt=frame,t=tdata)
    ax.set_title("Displacement versus time for $dt$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(3*np.pi,5*np.pi,100),interval = 100)
plt.show()