# Numerical Integration
#### Preliminaries
We have to import the array library `numpy` and the plotting library `matplotlib.pyplot`, note that we define shorter aliases for these.

Next we import from `numpy` some of the functions that we will use more frequently and from an utility library functions to format conveniently our results.

Eventually we tell to the plotting library how we want to format our plots (`fivethirtyeight` is a predefined style that mimics a popular site of statistical analysis and yes, statistical analysis and popular apply both to a single item... while `00_mplrc` refers to some tweaks I prefer to use)


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy import cos, exp, pi, sin, sqrt 
from IPython.display import Latex, display
plt.style.use(['fivethirtyeight', './00_mplrc'])
```

## Plot of the load

We start analyzing the loading process.

Substituting $t_1=1.2$s in the expression of $p(t)$ we have

$$p(t) = \left({15625 \over 27} {t^3 \over s^3} - {8750 \over 9} {t^2 \over s^2} 
      + {1100 \over 3} {t \over s} + 400 \right)\text{N}$$


```python
N = 6000
def it(t): return int((t+1)*N/3)

t = np.linspace(-1,2,N+1)
p = np.where(t<0, 0, np.where(
             t<=1.2, 15625./27*t**3 - 8750./9*t**2 + 1100./3*t + 400, 0))
```

We will use the above expression everywhere in the following.

Let's visualize this loading with the help of a plot:


```python
plt.plot(t,p)
plt.ylim((-20,520))
plt.xticks((-1,0,1,1.2,2))
plt.xlabel('t [second]')
plt.ylabel('p(t) [Newton]')
plt.title(r'The load $p(t)$')
plt.show()
```

From the plot it is apparent that the loading can be approximated by a constant rectangle plus a sine, of relatively small amplitude, with a period of about 1 second.

I expect that the forced response will be a damped oscillation about the static displacement $400/k$, followed by a free response in which we will observe a damped oscillation about the zero displacement.

## The Particular Integral

Writing the particuar integral as

$$\xi(t) = \left(A {t^3 \over s^3} + B {t^2 \over s^2} 
      + C {t \over s} + D \right)$$
      
where $A, B, C, D$ all have the dimension of a length, deriving with respect to time and substituting in the equation of motion we have

$$\left( 6A\frac{t}{s} + 2 B\right)\frac{m}{s^2} 
      + \left( 3A\frac{t^2}{s^2} + 2B\frac{t}{s} + C \right) \frac{c}{s}
      + \left( A\frac{t^3}{s^3} + B\frac{t^2}{s^2} + C\frac{t}{s} + D \right) k=\left({15625 \over 27} {t^3 \over s^3} - {8750 \over 9} {t^2 \over s^2} 
      + {1100 \over 3} {t \over s} + 400 \right)\text{N}$$
      
Collecting the powers of $t/s$ in the left member we have

$$A k \frac{t^3}{s^3} + \left( B k + \frac{3 A c}{s} \right) \frac{t^2}{s^2} +  \left( C k + \frac{2 B c}{s} + \frac{6 A m}{s^2} \right) \frac{t}{s}
      + \left(D k + \frac{C c}{s} + \frac{2 B m}{s^2}\right)=\left({15625 \over 27} {t^3 \over s^3} - {8750 \over 9} {t^2 \over s^2} 
      + {1100 \over 3} {t \over s} + 400 \right)\text{N}$$

and equating the coefficients of the different powers, starting from the higher to the lower ones, we have a sequence of equations of a single unknown, solving this sequence of equations gives, upon substitution of the numerical values of the system parameters,

\begin{align*}
  A& = \frac{625}{14256}\text{m}=0.0438411896745
\text{m}\\
  B& = \left(\frac{-175}{2376} - \frac{19}{209088} \sqrt{11}\right)\text{m}=-0.0739545830991
\text{m}\\
  C& = \left(\frac{133}{1306800} \sqrt{11} + \frac{218117}{7920000}\right)\text{m}=0.0278775758319\text{m}\\
  D& = \left(\frac{3013181}{99000000} - \frac{135571859}{7187400000000} \sqrt{11}\right)\text{m}=0.0303736121006\text{m}\\
\end{align*}

Substituting in $\xi(t)$ and taking the time derivative we have

\begin{align*}
  \xi(t) &= \frac{625}{14256} t^3 - 
            \left( \frac{19}{209088} \sqrt{11} + \frac{175}{2376} \right) t^2 +
            \left( \frac{133}{1306800} \sqrt{11} + \frac{218117}{7920000} \right) t -
            \frac{135571859}{7187400000000} \sqrt{11} + \frac{3013181}{99000000},\\
  \dot\xi(t) &= \frac{625}{4752} t^2 - 
              2 \left( \frac{19}{209088} \sqrt{11} + \frac{175}{2376} \right) t + 
              \frac{133}{1306800} \sqrt{11} + \frac{218117}{7920000}
\end{align*}


```python
xi = np.where(t<0, 0,
     np.where(t<=1.2001,
              625./14256*t**3 - (19./209088*sqrt(11) + 175./2376)*t**2 
          + t*(133./1306800*sqrt(11) + 218117./7920000) 
          - 135571859./7187400000000*sqrt(11) + 3013181./99000000, 0))
dot_xi = np.where(t<0, 0,
         np.where(t<=1.2001, 
                  625./4752*t**2 - 2*(t*(19./209088*sqrt(11) + 175./2376)) 
                + 133./1306800*sqrt(11) + 218117./7920000, 0))
xi_0 = - 135571859./7187400000000*sqrt(11) + 3013181./99000000
dot_xi_0 = + 133./1306800*sqrt(11) + 218117./7920000
```

### Plot of the particular integral


```python
plt.plot(t,xi)
plt.title('Particular integral')
plt.xlabel('s')
plt.ylabel('m')
plt.show()
plt.plot(t, dot_xi)
plt.title('Time derivative of particular integral')
plt.xlabel('s')
plt.ylabel('m/s')
plt.show()
```

## System Response

### Forced Response $0\le t\le t_1$

The response in terms of displacement and velocity is

\begin{align}
  x(t) &= \exp(-\zeta\omega_nt) \left(
    R\cos(\omega_Dt) + S\sin(\omega_Dt)\right) + \xi(t),\\
  v(t) &= \exp(-\zeta\omega_nt) \left(
    \left(S\cos(\omega_Dt)-R\sin(\omega_Dt)\right)\omega_D - 
    \left(R\cos(\omega_Dt)+S\sin(\omega_Dt)\right)\zeta\omega_n 
  \right) + \dot\xi(t)
\end{align}

and we can write the following initial conditions, taking into account that at time $t=0$ the system is at rest,

\begin{align}
  x(0) &= 1 \cdot \left(
    R\cdot1 + S\cdot0\right) + \xi_0\\&=R+\xi_0=0,\\
  \dot x(0) &= 1 \cdot \left(
    \left(S\cdot1-R\cdot0\right)\omega_D - 
    \left(R\cdot1+S\cdot0\right)\zeta\omega_n 
  \right) + \dot\xi_0\\&=S\omega_D-R\zeta\omega_n+\dot\xi_0=0
\end{align}

The constants of integration are

\begin{align}
  R &= -\xi_0,\\
  S &= \frac{R\zeta\omega_n-\dot\xi_0}{\omega_D}
\end{align}


```python
mass = 12.0
stif = 13200.0
zeta = 0.038

w2_n = stif/mass
w1_n = sqrt(w2_n)
w1_D = w1_n*sqrt(1-zeta**2)

damp = 2*zeta*w1_n*mass

R = -xi_0
S = (R*zeta*w1_n-dot_xi_0)/w1_D

display(HTML("<center><h3>Forced Displacements</h3></center>"))
display(Latex(r"""
    $x(t) = \exp(-%g\cdot%g\,t) \left(
    %+g\cos(%g\,t)%+g\sin(%g\,t)
    \right)+\xi(t)$
    """%(zeta,w1_n,R,w1_D,S,w1_D)))

t1  = 1.2
it1 = it(t1)

x_t1 = exp(-zeta*w1_n*t1)*(R*cos(w1_D*t1)+S*sin(w1_D*t1))+xi[it1]
v_t1 =(exp(-zeta*w1_n*t1)*(R*cos(w1_D*t1)+S*sin(w1_D*t1))*(-zeta)*w1_n
      +exp(-zeta*w1_n*t1)*(S*cos(w1_D*t1)-R*sin(w1_D*t1))*w1_D
      +dot_xi[it1])
display(Latex(
        r'$$x(t_1)=x_1=%+g\,\text{m},\qquad v(t_1)=v_1=%+g\,\text{m/s}$$'%
        (x_t1,v_t1)))
```


<center><h3>Forced Displacements</h3></center>




    $x(t) = \exp(-0.038\cdot33.1662\,t) \left(
    -0.0303736\cos(33.1423\,t)-0.00199618\sin(33.1423\,t)
    \right)+\xi(t)$
    



$$x(t_1)=x_1=+0.035918\,\text{m},\qquad v(t_1)=v_1=+0.237819\,\text{m/s}$$


### Free Response

For $t\ge t_1$ the response (no external force) is given by

\begin{align}
  x^*(t) &= \exp(-\zeta\omega_nt) \left(R^*\cos(\omega_Dt) + S^*\sin(\omega_Dt)\right),\\
  v^*(t) &= \exp(-\zeta\omega_nt) \left(
    \left(S^*\cos(\omega_Dt)-R^*\sin(\omega_Dt)\right)\omega_D - 
    \left(R^*\cos(\omega_Dt)+S^*\sin(\omega_Dt)\right)\zeta\omega_n 
  \right).
\end{align}

By the new initial conditions,

$$x^*(t_1) = x(t_1) = x_1, \qquad v^*(t_1) = v(t_1) = v_1,$$

we have, with 
$e_1 = \exp(-\zeta\omega_{n}t_1)$, 
$c_1 = \cos(\omega_Dt_1)$ and
$s_1 = \sin(\omega_Dt_1)$

$$ e_1 \begin{bmatrix} 
     c_1 & s_1 \\
     -\omega_D s_1 - \zeta\omega_n c_1 & \omega_D c_1 - \zeta\omega_n s_1 \end{bmatrix} \,
     \begin{Bmatrix}R^*\\S^*\end{Bmatrix} = \begin{Bmatrix}x_1\\v_1\end{Bmatrix}
$$

that gives

$$
\begin{Bmatrix}R^*\\S^*\end{Bmatrix} =
\frac{1}{\omega_D\,e_1}\,
\begin{bmatrix}
  \omega_D c_1 - \zeta\omega_n s_1 & -s_1 \\
  \zeta\omega_n c_1 + \omega_D s_1 & c_1 
\end{bmatrix}\,
\begin{Bmatrix}x_1\\v_1\end{Bmatrix}.
$$


```python
e_t1 = exp(-zeta*w1_n*t1)
c_t1 = cos(w1_D*t1)
s_t1 = sin(w1_D*t1)

Rs = (w1_D*c_t1*x_t1 - zeta*w1_n*s_t1*x_t1 - s_t1*v_t1) /w1_D /e_t1
Ss = (w1_D*s_t1*x_t1 + zeta*w1_n*c_t1*x_t1 + c_t1*v_t1) /w1_D /e_t1

display(HTML("<center><h3>Free Displacements</h3></center>"))
display(Latex(r"""
    $$x^*(t) = \exp(-%g\cdot%g\,t) \left(
    %+g\cos(%g\,t)%+g\sin(%g\,t)
    \right)$$
    """%(zeta,w1_n,Rs,w1_D,Ss,w1_D)))

xs_t1 = exp(-zeta*w1_n*t1)*(Rs*cos(w1_D*t1)+Ss*sin(w1_D*t1))
vs_t1 = ((exp(-zeta*w1_n*t1)*(Rs*cos(w1_D*t1)+Ss*sin(w1_D*t1))*(-zeta)*w1_n
      +exp(-zeta*w1_n*t1)*(Ss*cos(w1_D*t1)-Rs*sin(w1_D*t1))*w1_D))

display(Latex(r'$$x^*(t_1)=%+g\,\text{m},\qquad v^*(t_1) = %+g\,\text{m/s}$$'%
              (xs_t1,vs_t1)))
```


<center><h3>Free Displacements</h3></center>




    $$x^*(t) = \exp(-0.038\cdot33.1662\,t) \left(
    -0.112254\cos(33.1423\,t)+0.124351\sin(33.1423\,t)
    \right)$$
    



$$x^*(t_1)=+0.035918\,\text{m},\qquad v^*(t_1) = +0.237819\,\text{m/s}$$


#### Putting it all together

First, the homogeneous response


```python
x_hom = np.where(t<0, 0,
        np.where(t<1.2001, exp(-zeta*w1_n*t) * (R *cos(w1_D*t) + S *sin(w1_D*t)),
                           exp(-zeta*w1_n*t) * (Rs*cos(w1_D*t) + Ss*sin(w1_D*t))))
v_hom = np.where(t<0, 0,
        np.where(t<1.2001,
                 exp(-t*zeta*w1_n) * (
                    ( S*w1_D- R*zeta*w1_n)*cos(t*w1_D) -
                    ( S*zeta*w1_n+ R*w1_D)*sin(t*w1_D)),
                 exp(-t*zeta*w1_n) * (
                    (Ss*w1_D-Rs*zeta*w1_n)*cos(t*w1_D) -
                    (Ss*zeta*w1_n+Rs*w1_D)*sin(t*w1_D))))
```

then, we put together the homogeneous response and the particular integral, in the different intervals


```python
x = x_hom+xi
v = v_hom+dot_xi
```

### Plot of the response


```python
plt.plot(t,x)
plt.title('Displacement')
plt.show()

plt.plot(t,v)
plt.title('Velocity')
plt.show()

plt.plot(t,x)
plt.title('Displacement, zoom around t_1')
plt.xlim((1.16,1.24))
plt.show()

plt.plot(t,v)
plt.title('Velocity, zoom around t_1')
plt.xlim((1.16,1.24))
plt.show()
```

## Numerical Integration

The time step we are going to use is specified in terms of the natural period of vibration, $h=T_n/12$:


```python
t_n = 2*pi/w1_n
h = t_n/12.0
```

We need a function that returns the value of the load,


```python
def load(t):
    return 0 if t > 1.20001 else (
        + 578.703703703703705*t**3 
        - 972.22222222222222*t**2
        + 366.666666666666667*t
        + 400.)
```

### Initialization

The final time for the computation, the factors that modify the increment of the load to take into account the initial conditions, the modified stiffness, the containers for the results, the initial values of the results.


```python
stop = 2.0 + h/2

a0fac = 3.0*mass + damp*h/2.0
v0fac = 6.0*mass/h + 3.0*damp

ks = stif + 3.0*damp/h + 6.0*mass/h**2

T, X, V, A = [], [], [], []

x0 = 0
v0 = 0
t0 = 0
p0 = load(t0)
```

### Computational Loop


```python
while t0 < stop:
    a0 = (p0-stif*x0-damp*v0)/mass
    for current, vec in zip((t0,x0,v0,a0),(T,X,V,A)):
        vec.append(current)
    t1 = t0+h
    p1 = load(t1)
    dp = p1-p0
    dx = (dp+a0fac*a0+v0fac*v0)/ks
    dv = 3*dx/h -a0*h/2 -3*v0
    t0 = t1 ; p0 = p1 ; x0 = x0+dx ; v0 = v0+dv
```

### The Plot


```python
plt.plot(T,X,label='numerical')
plt.plot(t,x,label='analytical')
plt.xlim((0,2))
plt.legend(loc=3);
```


```python
plt.plot(T,X,'-o',label='numerical')
plt.plot(t,x,'-x',label='analytical')
plt.title('Displacement, zoom around t_1')
plt.xlim((1.16,1.24))
plt.ylim((0,0.04))
plt.legend(loc=3)
plt.show()
```


```python
# an IPython incantation that properly formats this notebook
from IPython.display import HTML
HTML(open("00_custom.sav.css", "r").read())
```




<style>
@font-face {
    font-family: 'Charis SIL';
    src: url('fonts/CharisSILEur-R.woff') format('woff');
}
@font-face {
    font-family: 'Architects Daughter';
    font-style: normal;
    src: url(https://fonts.gstatic.com/s/architectsdaughter/v6/RXTgOOQ9AAtaVOHxx0IUBM3t7GjCYufj5TXV5VnA2p8.woff2) format('woff2');
    unicode-range: U+0000-00FF, U+0131, U+0152-0153, U+02C6, U+02DA, U+02DC, U+2000-206F, U+2074, U+20AC, U+2212, U+2215, U+E0FF, U+EFFD, U+F000;
}
.prompt{display: None;}
div.cell{margin: auto;width:900px;}

div #notebook { /* centre the content */
    background: #fffaf0;
    margin: auto;
    padding: 0em;
    padding-top: 1em;
}
div #notebook_container {width: 960px!important;}
#notebook li { /* More space between bullet points */
    margin-top:0.2em;
}

div.text_cell_render{
    font-family: "Charis SIL", Cambria,  serif;
    line-height: 155%;
    font-size: 130%;
}
.CodeMirror{
    font-family: Consolas, monospace;
    font-size: 140%;
    background-color:#fcffff!important;
}

.output_area {
    font-family: Consolas,monospace;font-size: 120%;
    background-color:#fcffff!important;
    margin-top:0.8em;
}

div.output_latex{overflow:hidden}

h1 {
    font-family: 'Architects Daughter', serif;
    font-size:   48pt!important;
    text-align:  center;
    text-shadow: 4px 4px 4px #aaa;
    padding-bottom: 48pt;
}

.warning{color: rgb( 240, 20, 20 )}
</style>





cff4c799de5df2b5a1c2d3fdbed691ac91776afd
ipynb
Jupyter Notebook
dati_2015/ha03/03_Numerical_Integration.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2015/ha03/03_Numerical_Integration.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2015/ha03/03_Numerical_Integration.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
2019-06-23T12:32:39.000Z
2021-08-15T18:33:55.000Z
Qwen/Qwen-72B
1. YES
2. YES