# Solution {-}

Evaluation of integrals analytically can be quite laborious (some also impossible), so one should not overlook the possibility of using numerical integration. To demonstrate the effectiveness of the numerical approach, consider the following system driven by white noise whose PSD = 10 units.

\begin{equation*}
  G(s)=\frac{1}{s^2+\omega_0^2}, \quad \omega_0=20\pi
\end{equation*}

Lets say we want to find $E[x^2(t)]$ for $0 \leq t \leq 1$ with a sample spacing of 0.001 seconds.


```python
from numpy import arange, pi, sin, cos
from scipy.integrate import quad
from sympy import inverse_laplace_transform, integrate, symbols
import matplotlib.pyplot as plt

s, t, w0, A = symbols('s t w0 A', positive = True)

# Filter in s-domain
G = 1/(s**2 + w0**2)
display(G)
```


$\displaystyle \frac{1}{s^{2} + w_{0}^{2}}$



```python
# Filter in t-domain
g = inverse_laplace_transform(G, s, t)
display(g)
```


$\displaystyle \frac{\sin{\left(t w_{0} \right)}}{w_{0}}$



```python
EX2 = integrate(A*g**2, t)
display(EX2)
```


$\displaystyle \frac{A \left(\frac{t w_{0}}{2} - \frac{\sin{\left(t w_{0} \right)} \cos{\left(t w_{0} \right)}}{2}\right)}{w_{0}^{3}}$



```python
# Numerical values
dt = .001
A = 10
w0 = 20*pi

# Initialize plot vectors
exact = []
approx = []
time=arange(0, 0.1, dt)

# Exact integral
for i in range(100):
    EX2 = A/w0**3*(1/2*w0*i*dt - 1/2*sin(w0*i*dt)*cos(w0*i*dt))
    exact.append(EX2) 

# Numerical integration
f = lambda dt: A*(sin(w0*dt)/w0)**2
for i in range(100):
    EX2 = quad(f, 0, i*dt)
    approx.append(EX2[0])
    
plt.title('Square Impulsive Respons')
plt.xlabel('Time (second)')
plt.ylabel('E(X2)')
plt.plot(time, exact, 'b', label='Exact')
plt.plot(time, approx, 'r', label='Approx')
plt.legend(loc='lower right')
plt.grid(True, which='both')
plt.show()
```

33e6a056315978baa42ba9786c19b47f5a3160ad
ipynb
Jupyter Notebook
Problem 3.16.ipynb
mfkiwl/GMPE340
3602b8ba859a2c7db2cab96862472597dc1ac793
MIT
2021-09-07T09:36:36.000Z
2021-09-07T09:36:36.000Z
Problem 3.16.ipynb
mfkiwl/GMPE340
3602b8ba859a2c7db2cab96862472597dc1ac793
MIT
Problem 3.16.ipynb
mfkiwl/GMPE340
3602b8ba859a2c7db2cab96862472597dc1ac793
MIT
2021-09-20T18:48:20.000Z
2021-09-20T18:48:20.000Z
Qwen/Qwen-72B
1. YES
2. YES