#  Arbitrary Periodic Forcing (Fourier Series)

## Sympy

SymPy is a Python package for symbolic computing. It can do many symbolic operations, for instance, integration, differentiation, linear algebra, etc. See http://sympy.org for more details of the features and the documentation. Today we will cover how to do integrals using SymPy and use it to find the Fourier series that represents a sawtooth function.


```python
import sympy as sm
```

The function `init_printing()` enables LaTeX based rendering in the Jupyter notebook of all SymPy objects.


```python
sm.init_printing()
```

Symbols can be created by using the `symbols()` function.


```python
x, y, z = sm.symbols('x, y, z')
```


```python
x, y, z
```

The `integrate()` function allows you to do symbolic indefinite or definite integrals. Note that the constants of integration are not included in indefinite integrals.


```python
sm.integrate(x * y, x)
```

The `Integral` class creates and unevaluated integral, where as the `integrate()` function automatically evaluates the integral.


```python
expr = sm.Integral(x * y, x)
expr
```

To evaluate the unevaluated form you call the `.doit()` method. Note that all unevaluated SymPy objects have this method.


```python
expr.doit()
```

This shows how to create an unevaluated definite integral, store it in a variable, and then evaluate it.


```python
expr = sm.Integral(x * y, (x, 0, 5))
expr
```


```python
expr.doit()
```

# Fourier Coefficients for the Sawtooth function

Now let's compute the Fourier coefficients for a saw tooth function. The function that describes the saw tooth is:

$$
F(t) = 
\begin{cases} 
      A \left( \frac{4t}{T} - 1 \right) & 0 \leq t \leq T/2 \\
      A \left( 3 - \frac{4t}{t} \right) & T/2 \leq t \leq T 
\end{cases}
$$

where:

- $A$ is the amplitude of the saw tooth
- $T$ is the period of the saw tooth
- $\omega_T$ is the frequency of the saw tooth, i.e. $\omega_T = \frac{2\pi}{T}$
- $t$ is time

This is a piecewise function with two parts from $t=0$ to $t=T$.


```python
A, T, wT, t = sm.symbols('A, T, omega_T, t', real=True, positive=True)
A, T, wT, t
```

The first Fourier coefficient $a_0$ describes the average value of the periodic function. and is:

$$a_0 = \frac{2}{T} \int_0^T F(t) dt$$

This integral will have to be done in two parts:

$$a_0 = a_{01} + a_{02} = \frac{2}{T} \int_0^{T/2} F(t) dt + \frac{2}{T} \int_{T/2}^T F(t) dt$$

These two integrals are evaluated below. Note that $a_0$ evaluates to zero. This is because the average of our function is 0.


```python
ao_1 = 2 / T * sm.Integral(A * (4 * t / T - 1), (t, 0, T / 2))
ao_1
```


```python
ao_1.doit()
```


```python
ao_2 = 2 / T * sm.Integral(A * (3 - 4 * t / T), (t, T / 2, T))
ao_2
```


```python
ao_2.doit()
```

But SymPy can also handle piecewise directly. The following shows how to define a piecewise function.


```python
F_1 = A * (4 * t / T - 1)
F_2 = A * (3 - 4 * t / T)
```


```python
F = sm.Piecewise((F_1, t<=T/2),
                 (F_2, T/2<=t))
F
```

The integral can be taken of the entire piecewise function in one call.


```python
sm.integrate(F, (t, 0, T))
```

Now the Fourier coefficients $a_n$ and $b_n$ can be computed.

$$
a_n = \frac{2}{T}\int_0^T F(t) \cos n\omega_Tt dt \\
b_n = \frac{2}{T}\int_0^T F(t) \sin n\omega_Tt dt
$$


```python
n = sm.symbols('n', real=True, positive=True)
```

For $a_n$:


```python
an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))
an
```


```python
an.doit()
```

This can be simplified:


```python
an = an.doit().simplify()
an
```

Now substitute the $2\pi/T$ for $\omega_T$.


```python
an = an.subs({wT: 2 * sm.pi / T})
an
```

Let's see how this function varies with increasing $n$. We will use a loop but the SymPy expressions will not automatically display because they are inside a loop. So we need to use SymPy's `latex()` function and the IPython display tools. SymPy's `latex()` function  transforms the SymPy expression into a string of matching LaTeX commands.


```python
sm.latex(an, mode='inline')
```




    '$\\frac{A}{2 \\pi^{2} n^{2}} \\left(- 2 \\pi n \\sin{\\left (2 \\pi n \\right )} + 8 \\cos{\\left (\\pi n \\right )} - 4 \\cos{\\left (2 \\pi n \\right )} - 4\\right)$'



The `display()` and `LaTeX()` functions then turn the LaTeX string in to a displayed version.


```python
from IPython.display import display, Latex 
```

Now we can see how $a_n$ varies with $n=1,2,\ldots$.


```python
for n_i in range(1, 6):
    ans = an.subs({n: n_i})
    display(Latex('$a_{} = $'.format(n_i) + sm.latex(ans, mode='inline')))
```


$a_1 = $$- \frac{8 A}{\pi^{2}}$



$a_2 = $$0$



$a_3 = $$- \frac{8 A}{9 \pi^{2}}$



$a_4 = $$0$



$a_5 = $$- \frac{8 A}{25 \pi^{2}}$


For even $n$ values the coefficient is zero and for even values it variers with the inverse of $n^2$. More precisely:

$$
a_n =
\begin{cases}
0 & \textrm{if }n\textrm{ is even} \\
-\frac{8A}{n^2\pi^2} & \textrm{if }n\textrm{ is odd}
\end{cases}
$$

SymPy can actually reduce this further if your set the assumption that $n$ is an integer.


```python
n = sm.symbols('n', real=True, positive=True, integer=True)
an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))
an = an.doit().simplify()
an.subs({wT: 2 * sm.pi / T})
```

The odd and even versions can be computed by setting the respective assumptions.


```python
n = sm.symbols('n', real=True, positive=True, integer=True, odd=True)
an = 2 / T * sm.Integral(F * sm.cos(n * wT * t), (t, 0, T))
an = an.doit().simplify()
an.subs({wT: 2 * sm.pi / T})
```

Note that $b_n$ is always zero:


```python
bn = 2 / T * sm.Integral(F * sm.sin(n * wT * t), (t, 0, T))
bn
```


```python
bn.doit().simplify().subs({wT: 2 * sm.pi / T})
```

# Random questions questions from class

The source code for any function is shown with the double question mark.


```python
sm.Integral??
```

The location of the file the function comes from can be found by accessing the `__file__` attribute.


```python
sm.__file__
```




    '/home/moorepants/miniconda3/lib/python3.5/site-packages/sympy/__init__.py'



SymPy can handle unit for you too.


```python
from sympy.physics.units import meter, kilogram, second
```


```python
a = 10 * second
a
```


```python
b = 5 * meter
```


```python
c = 6 * kilogram
```


```python
a, b, c
```


```python
expr = c * b / a**2
expr
```


```python
from sympy.physics.units import newton
```


```python
5 * newton
```

# Numerical evalution of the Fourier Series

Now the Fourier coefficients can be used to plot the approximation of the saw tooth forcing function.


```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
```

The following function plots the actual sawtooth function. It does it all in one line by cleverly using the absolute value and the modulo functions.


```python
def sawtooth(A, T, t):
    return (4 * A / T) * (T / 2 - np.abs(t % T - T / 2) ) - A
```


```python
A = 1
T = 2
t = np.linspace(0, 5, num=500)
```


```python
plt.figure()

plt.plot(t, sawtooth(A, T, t));
```


    <IPython.core.display.Javascript object>






# Exercise

Write a function that computes the Fourier approximation of the sawtooth function for a given value of $n$, i.e. using a finite number of terms. Then plot it for $n=2, 4, 6, 8, 10$ on top of the actual sawtooth function. How many terms of the infinite series are needed to get a good sawtooth?


```python
def sawtooth_approximation(n, A, T, t):
    # type code here
    return    
```

# Exercise

Now, using the Fourier series you've developed compute the response of a linear system to the sawtooth function and plot the steady state trajectory of $x$ give:

$$m\ddot{x} + c\dot{x} + kx = F(t)$$

Use:

- $m$: 100 kg
- $\omega_n$: 10 rad/s
- $\zeta$: 0.01
- $T=\frac{2\pi}{5 rad/s}$
- $A=1000N$

Investigate the difference in steady state response for different approximations of $F(t)$, i.e. various $n$, and report your findings in a markdown cell.


```python

```

5c1cc1b8ec58b9998666543e7ff9fd223b2c3f4d
ipynb
Jupyter Notebook
content/materials/notebooks/sawtooth_forcing.ipynb
greg1877/eng122
3a9821f41c860e190cd3c8285fac212ae897337b
CC-BY-4.0
2020-05-11T05:56:54.000Z
2020-12-20T13:44:12.000Z
content/materials/notebooks/sawtooth_forcing.ipynb
greg1877/eng122
3a9821f41c860e190cd3c8285fac212ae897337b
CC-BY-4.0
2016-09-20T01:26:33.000Z
2020-01-22T02:10:33.000Z
content/materials/notebooks/sawtooth_forcing.ipynb
greg1877/eng122
3a9821f41c860e190cd3c8285fac212ae897337b
CC-BY-4.0
2017-06-09T01:28:45.000Z
2020-02-29T02:36:18.000Z
Qwen/Qwen-72B
1. YES
2. YES