This is an IPyton notebook and it is coded in the Python programming language.

Python has the idea of modules, separate units of code that provide the programmer with ready to use, tested code. You can write yourself your modules, or rather use library modules, in any case you access functions or data defined in modules import-ing them.

For example, in the next code cell we import three functions from the math module (and a constant too) and then we do some other stuff...

import functions and data from math,
import all the stuff from sympy (symbolic algebra) under the prefix sy, and use this convention to initialize the sympy printing system,
import from the submodule IPython.display two utility functions,
eventually import the module style and execute its function clean that changes the default appearance of a notebook using some css styling and a small sprinkle of javascript.



In [1]:

    
%pylab --no-import-all inline

from math import exp, log, pi, sqrt

import sympy as sy ; sy.init_printing(use_latex=True)
from IPython.display import display, Latex

import style ; style.clean()









    



Populating the interactive namespace from numpy and matplotlib






    Out[1]:

Damped response for a given $x_0$

What is the analytical expression of the response of a damped system subject to the (quite common) initial conditions $x(0)=x_0$ and $\dot x(0)=0$?

We start by creating a set of symbols and giving them good names...



In [2]:

    
A, B, t, z, wn, wd, x0 = sy.symbols(
    'A, B, t, zeta, omega_n, omega_d, x_0')
display((A, B, t, z, wn, wd, x0))









    





$$\begin{pmatrix}A, & B, & t, & \zeta, & \omega_{n}, & \omega_{d}, & x_{0}\end{pmatrix}$$

Then we give the name u to the symbolic representation of the displacement, and the name v to the first time derivative of u.



In [3]:

    
u = sy.exp(-z*wn*t)*(A*sy.cos(wd*t)+B*sy.sin(wd*t))
v = u.diff(t).expand().collect(sy.exp(-z*wn*t))
v = v.collect(sy.cos(wd*t)).collect(sy.sin(wd*t))
display(Latex(r'\begin{align*}x(t)&='+sy.latex(u)+r'''\\
\\dot{x}(t) &= '''+sy.latex(v)+r'\end{align*}'))









    





\begin{align*}x(t)&=\left(A \cos{\left (\omega_{d} t \right )} + B \sin{\left (\omega_{d} t \right )}\right) e^{- \omega_{n} t \zeta}\\\dot{x}(t) &= \left(\left(- A \omega_{d} - B \omega_{n} \zeta\right) \sin{\left (\omega_{d} t \right )} + \left(- A \omega_{n} \zeta + B \omega_{d}\right) \cos{\left (\omega_{d} t \right )}\right) e^{- \omega_{n} t \zeta}\end{align*}

Now,

we form an equation substituting $t=0$ in u and v and subtracting $x_0$ from the first equation
we ask our sympy aka sy library to solve the system of equations



In [4]:

    
eq = (u.subs(t,0)-x0, v.subs(t,0))
coeff = sy.solve(eq, A,B)
display(coeff)
display(u.subs(coeff).collect(x0))
display(u.subs(coeff).collect(x0).subs(wn,wd/sy.sqrt(1-z*z)))









    





$$\begin{Bmatrix}A : x_{0}, & B : \frac{\omega_{n} x_{0}}{\omega_{d}} \zeta\end{Bmatrix}$$






    





$$x_{0} \left(\cos{\left (\omega_{d} t \right )} + \frac{\omega_{n} \zeta}{\omega_{d}} \sin{\left (\omega_{d} t \right )}\right) e^{- \omega_{n} t \zeta}$$






    





$$x_{0} \left(\frac{\zeta \sin{\left (\omega_{d} t \right )}}{\sqrt{- \zeta^{2} + 1}} + \cos{\left (\omega_{d} t \right )}\right) e^{- \frac{\omega_{d} t \zeta}{\sqrt{- \zeta^{2} + 1}}}$$

To plot the response we need to get rid of all the symbols in our expression, except for $t$, so we choose an arbitrary set of values for the free variables...



In [5]:

    
values = {x0:1.0, z:0.05, wd:2*sy.pi, wn:2*sy.pi/sy.sqrt(1-z*z)}
values[sy.pi] = sy.pi.evalf()
display(values)









    





$$\begin{Bmatrix}\pi : 3.14159265358979, & \omega_{d} : 2 \pi, & \omega_{n} : \frac{2 \pi}{\sqrt{- \zeta^{2} + 1}}, & x_{0} : 1.0, & \zeta : 0.05\end{Bmatrix}$$

At the end, we display the analytical expression of $x(t)$ and plot

$x(t)$,
$\dot x(t)$.



In [6]:

    
display(u.subs(coeff).subs(values))
sy.plot(u.subs(coeff).subs(values), (t,0,5.99));
sy.plot(v.subs(coeff).subs(values), (t,0,5.99)  );









    





$$\left(0.0500626174321759 \sin{\left (2 \pi t \right )} + 1.0 \cos{\left (2 \pi t \right )}\right) e^{- 0.100125234864352 \pi t}$$

Logarithmic decrement

We assume that our system has been deformed by the application of some static load, so that at time $t=0$ the initial conditions are $x(0)=x_0$ and $\dot x(0)=0$.

At $t=0$ the static force is suddenly released, we observe an oscillatory behaviour as we have detailed in the preceding section.

and we are able to measure the peak displacement after $n$ cycles, $x_n$.

Specify data

The following code gives appropriate names to values measured or computed from expressions, the value of the logarithmic decrement is eventually printed.



In [7]:

    
n = 3
x0 = 8.00
xn = 3.11
dn = log(x0/xn)
print "delta_n =",dn









    



delta_n = 0.944818815489

The exact expression for $\zeta$ is

$$ \zeta = \sqrt{\frac{\delta_n^2}{(2 n \pi)^2 + \delta_n^2}} \text{ with } \delta_n = \log\frac{x_0}{x_n}.$$

We've already computed $\delta_n$ in the previous code cell so we can write the numerator and the denominator of the radicand and eventually compute the "exact" value of $\zeta$.



In [8]:

    
num = dn*dn
den = 4 * n*n * pi*pi + dn*dn
z = sqrt(num/den)
print "zeta = %f"%(z*100,)









    



zeta = 5.006135

Approximate expression

The linearized expression for $\zeta$, that we have found using the first and second terms of the McLaurin series expressing $\exp(a\zeta)=1+a\zeta$, is

$$\zeta \approx \frac{1}{2 n \pi}\,\frac{x_0-x_n}{x_n}.$$

We know that this approximation is good for low damping and few cycles...



In [9]:

    
print "Estimate", 100 * (x0-xn) / xn / ( 2*n* pi )
print "error (z/est.)",z /((x0-xn) / xn / ( 2*n* pi ))









    



Estimate 8.34156132604
error (z/est.) 0.600143600147

Iterative solution

We write, formally,

$$\zeta = \xi\sqrt{1-\zeta^2}, \qquad\text{with }\xi = \frac{1}{2 n \pi} \log\frac{x_0}{x_n}$$

and derive an iterative formulation

$$\zeta_{i+1} = \xi\sqrt{1-\zeta_i^2}, \qquad\text{with } i = 0, 1, 2, \ldots$$

Note that, for $\zeta_0 = 0$, it is $\zeta_1 = \xi$.

In theory, we have to check the convergence of the sequence and stop when the convergence is reached, but in our cases the convergence is soo fast that we are safe implementing the solution in a for loop



In [10]:

    
csi= dn/2/n/pi
z_i = csi

print "      zeta_0 = 0.00000000000000%"
print "csi = zeta_1 = %16.14f%%" % (csi*100,)
for i in range(2, 6):
    z_i = csi*sqrt(1-z_i*z_i)
    print '      zeta_%d = %16.14f%%'%(i,100*z_i)









    



      zeta_0 = 0.00000000000000%
csi = zeta_1 = 5.01241949370850%
      zeta_2 = 5.00611884484703%
      zeta_3 = 5.00613468475389%
      zeta_4 = 5.00613464495722%
      zeta_5 = 5.00613464505721%