This notebook is coded in python, a programmming language built around the idea of modules where most of the code resides.

You can write your modules or use library modules, in any case you access functions or classes defined in modules, import-ing them from the modules.

That is said to explain why we import three functions from the math module (and a constant too), and some other stuff too...



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

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)
display(u)









    





$$\left(A \cos{\left (\omega_{d} t \right )} + B \sin{\left (\omega_{d} t \right )}\right) e^{- \omega_{n} t \zeta}$$

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 [5]:

    
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 [6]:

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









    





$$\begin{Bmatrix}\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 [7]:

    
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 [8]:

    
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 [9]:

    
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 [10]:

    
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}, \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}, \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 [11]:

    
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%



In [ ]: