Taylor Examples 5.4 and 5.5
Contents
Taylor Examples 5.4 and 5.5#
last revised: 17-Jan-2019 by Dick Furnstahl [furnstahl.1@osu.edu]
Here we are exploring the Fourier series for a waveform defined to be even about the origin, so \(f(-t) = +f(t)\), with period \(\tau\). That means that the integrand for the \(b_m\) coefficients is odd and so all of the corresponding integrals vanish.
The particular wave of interest here is a periodic rectangular pulse, such that in the interval \(-\tau/2 \leq t \leq \tau/2\), the function takes the form:
\(\newcommand{\fmax}{f_{\textrm{max}}}\) \(\begin{align} f(t) = \left\{ \begin{array}{lc} \fmax & -\Delta\tau/2 < t < \Delta\tau/2 \\ 0 & |t| > \Delta\tau/2 \end{array} \right. \end{align}\)
As usual, we do the \(a_0\) integral separately and then \(a_m\) for \(m \geq 1\). Here \(\omega \equiv 2\pi/\tau\). The results are (also see the text):
\(\begin{align} a_0 = \frac{1}{\tau} \int_{-\tau/2}^{\tau/2} f(t)\, dt = \frac{\fmax\Delta\tau}{\tau} \end{align}\)
\(\begin{align} a_m = \frac{2}{\tau} \int_{-\tau/2}^{\tau/2} \cos(m\omega t) f(t)\, dt = \frac{4\fmax}{\tau} \int_0^{\Delta\tau/2} \cos(2\pi m t/\tau) \, dt % &= - \frac{2\fmax}{m\pi}\left[\cos(m\pi t)\right]^1_0 \\ = \frac{2 \fmax}{m\pi} \sin\left(\frac{\pi m \Delta\tau}{\tau}\right) \end{align}\)
Define the functions we’ll need#
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
We start by defining a function for the rectangular pulse at any \(t\). The definition here is for a scalar function. That is, it won’t work to call it with \(t\) and array of time points, unlike other functions we have defined. It is possible to make it work, but then the function will be much less clear. When we need to evaluate it for all elements of an array t_pts, we will use the construction:
np.array([rectangular(t) for t in t_pts]).
def rectangular(t, tau=1., Delta_tau=0.25, f_max=1.):
"""Returns the triangular wave of amplitude f_max and even about the
origin at time t. The period is tau. It is defined as a scalar
function (i.e., only one value of t can be passed at a time).
**There must be a better way to define this!**
"""
t_adjust = t - round(t) # shift t so we can test for being within Delta_tau
if np.fabs(t_adjust) < Delta_tau/2.:
return f_max
else:
return 0.
Now a function that creates an array of Fourier coefficients for the triangular wave up to order N_max.
def rectangular_coeffs_by_hand(N_max, tau=1., Delta_tau=0.25, f_max=1):
"""Fourier coefficients calculated by hand and loaded into an array
as pairs in the form (a_n, b_n).
"""
# initialize the array with the vaues of a_0 and b_0
coeffs_array = [(f_max * Delta_tau / tau, 0)]
for n in np.arange(1, N_max, 1):
# these are pre-calcualated values of a_n and b_n
a_n = 2.*f_max/(n * np.pi) * np.sin((np.pi * n * Delta_tau) / tau)
b_n = 0.
coeffs_array.append((a_n, b_n)) # add them to the end of the array
return np.array(coeffs_array) # convert to a numpy array to be safe
We would like a general way to construct the array of Fourier coefficients given any periodic function. Our first pass at that uses a class definition and the scipy integration function quad.
class FourierSeries():
"""
Fourier series class finds the coefficients in a Fourier series with
period tau up to a specified order.
Assume these imports:
from scipy.integrate import quad
import numpy as np
***Need to add args to the function.***
"""
def __init__(self,
function,
tau=1.,
N_max=10
):
self.function = function
self.tau = tau
self.omega = 2. * np.pi / tau
self.N_max = N_max
# add something to quit if Nmax < 0 or not an integer (try and except)
def a0_calc(self):
"""Calculate the constant Fourier coefficient by integration"""
answer, error = quad(self.function, -tau/2., tau/2., args=(tau,))
return (1./self.tau) * answer
def an_integrand(self, t, n):
"""Integrand for the nth cosine coefficient"""
return self.function(t, tau) * np.cos(n * self.omega * t)
def an_calc(self, n):
"""Calculate the nth cosine coefficient (n > 0)"""
# note comma after n in args
answer, error = quad(self.an_integrand, -tau/2., tau/2., args=(n,))
return (2./self.tau) * answer
def bn_integrand(self, t, n):
"""Integrand for the nth cosine coefficient"""
return self.function(t, tau) * np.sin(n * self.omega * t)
def bn_calc(self, n):
"""Calculate the nth cosine coefficient (n > 0)"""
answer, error = quad(self.bn_integrand, -tau/2., tau/2., args=(n,))
return (2./self.tau) * answer
def coeffs_upto_Nmax(self):
"""Calculate the Fourier series up to Nmax"""
# first generate the coefficient
coeffs_array = [(self.a0_calc(), 0)] # a_0 and b_0
for n in np.arange(1, N_max, 1):
a_n = self.an_calc(n)
b_n = self.bn_calc(n)
coeffs_array.append((a_n, b_n)) # append a tuple of coefficients
return np.array(coeffs_array) # convert to a numpy array
Finally, we need a function that can take as input an array of t values and an array of Fourier coefficients and return the function at those t values with terms up to order N_max.
def Fourier_reconstruct(t_pts, coeffs_array, tau, N_max):
"""Sum up the Fourier series up to n = N_max terms."""
omega = 2. * np.pi / tau # definition of omega
result = 0.
# iterate over coefficients but only up to N_max
for n, (a,b) in enumerate(coeffs_array[:N_max+1]):
result = result + a * np.cos(n * omega * t_pts) \
+ b * np.sin(n * omega * t_pts)
return result
Example 5.4#
Ok, now we can do Example 5.4. Calculate the coefficients both ways.
N_max = 20
tau = 1.
f_max = 1.
Delta_tau = 0.25
coeffs_by_hand = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
fs = FourierSeries(rectangular, tau, N_max)
coeffs_by_quad = fs.coeffs_upto_Nmax()
Let’s check that the exact and numerical calculation of the coefficients agree.
(Note the space in the formats below, e.g., {a1: .6f}. This means to leave an extra space for a positive number so that it aligns at the decimal point with negative numbers.)
print(' n a_exact a_quad b_exact b_quad')
for n, ((a1,b1), (a2,b2)) in enumerate(zip(coeffs_by_hand, coeffs_by_quad)):
print(f'{n:2d} {a1: .6f} {a2: .6f} {b1: .6f} {b2: .6f}')
n a_exact a_quad b_exact b_quad
0 0.250000 0.250000 0.000000 0.000000
1 0.450158 0.450158 0.000000 0.000000
2 0.318310 0.318310 0.000000 0.000000
3 0.150053 0.150053 0.000000 0.000000
4 0.000000 0.000000 0.000000 0.000000
5 -0.090032 -0.090032 0.000000 0.000000
6 -0.106103 -0.106103 0.000000 0.000000
7 -0.064308 -0.064308 0.000000 0.000000
8 -0.000000 0.000000 0.000000 0.000000
9 0.050018 0.050018 0.000000 0.000000
10 0.063662 0.063662 0.000000 0.000000
11 0.040923 0.040923 0.000000 0.000000
12 0.000000 0.000000 0.000000 0.000000
13 -0.034628 -0.034628 0.000000 0.000000
14 -0.045473 -0.045473 0.000000 0.000000
15 -0.030011 -0.030011 0.000000 0.000000
16 -0.000000 0.000000 0.000000 0.000000
17 0.026480 0.026480 0.000000 0.000000
18 0.035368 0.035368 0.000000 0.000000
19 0.023693 0.023693 0.000000 0.000000
Make the comparison plot requested: N_max = 2 vs. N_max = 6.
t_pts = np.arange(-2., 4., .01)
f_pts_2 = Fourier_reconstruct(t_pts, coeffs_by_hand, tau, 2)
f_pts_6 = Fourier_reconstruct(t_pts, coeffs_by_hand, tau, 6)
# Python way to evaluate the rectangular function at an array of points:
# * np.array creates a numpy array;
# * note the []s around the inner statement;
# * rectangular(t) for t in t_pts
# means step through each element of t_pts, call it t, and
# evaluate rectangular at that t.
# * This is called a list comprehension. There are more compact ways,
# but this is clear and easy to debug.
rectangular_t_pts = np.array([rectangular(t, tau, Delta_tau, f_max) \
for t in t_pts])
# It looks like we should define a function for the axis!
fig_1 = plt.figure(figsize=(10,5))
ax_1 = fig_1.add_subplot(1,2,1)
ax_1.plot(t_pts, f_pts_2, label='N = 2', color='blue')
ax_1.plot(t_pts, rectangular_t_pts, label='exact', color='red')
ax_1.set_xlim(-1.1,4.1)
ax_1.set_xlabel('t')
ax_1.set_ylabel('f(t)')
ax_1.set_title('N = 2')
ax_1.legend()
ax_2 = fig_1.add_subplot(1,2,2)
ax_2.plot(t_pts, f_pts_6, label='N = 6', color='blue')
ax_2.plot(t_pts, rectangular_t_pts, label='exact', color='red')
ax_2.set_xlim(-1.1,4.1)
ax_2.set_xlabel('t')
ax_2.set_ylabel('f(t)')
ax_2.set_title('N = 6')
ax_2.legend();
fig_1.tight_layout()
fig_1.savefig('example_5.4.png')
Example 5.5#
We consider a driven oscillator with an even driving force so only cosine terms are non-zero in the Fourier transform of \(f(t)\). Linearity of the equation then means that the solution takes the form
\(\begin{align} x(t) = \sum_{n=0}^{\infty} A_n \cos(n\omega t - \delta_n) \;, \end{align}\)
where
\(\begin{align} A_n &= \frac{f_n}{\sqrt{(\omega_0^2 - n^2 \omega^2)^2 + 4\beta^2 n^2 \omega^2}} \\ \delta_n &= \arctan\left(\frac{2\beta n\omega}{\omega_0^2 - n^2 \omega^2} \right) \end{align}\)
For Example 5.5 we apply this to the rectangular driving force, so we can use the coefficients from Example 5.4, which play the role of the \(f_n\) coefficients here.
def driven_coefficients(tau_0, tau, beta, coeffs):
"""Calculate the coefficients A_n and delta_n for a driven oscillator."""
omega_0 = 2. * np.pi / tau_0
omega = 2. * np.pi / tau
A_n = []
delta_n = []
# the second entry in the tuple is called dummy because it is not used
for n, (f_n, dummy) in enumerate(coeffs):
A_n.append(f_n / np.sqrt((omega_0**2 - n**2 * omega**2)**2 +
(2. * beta * n *omega)**2))
delta_n.append(np.arctan2(2. * beta * n *omega,
omega_0**2 - n**2 * omega**2 ))
return np.array(A_n), np.array(delta_n)
def x_drive(t_pts, tau_0, tau, beta, N_max, coeffs):
"""Sum up the Fourier series for x(t) up to n = N_max terms."""
omega = 2. * np.pi / tau
A_n_array, delta_n_array = driven_coefficients(tau_0, tau, beta, coeffs)
result = 0.
# iterate over coefficients but only up to N_max
for n, (A_n, delta_n) in enumerate(zip(A_n_array[:N_max+1],
delta_n_array[:N_max+1])):
result = result + A_n * np.cos( n * omega * t_pts - delta_n)
return result
Ok, let’s do Example 5.5!
N_max = 20
tau = 1.
f_max = 1.
Delta_tau = 0.25
coeffs_by_hand = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
t_pts = np.arange(0., 7., .01)
# first set of parameters
N_max = 20
f_max = 1.
Delta_tau = 0.25
tau = 1.
coeffs_by_hand_1 = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
tau_0_1 = 1.
beta = 0.2
A_n_array_1, delta_n_array_1 = driven_coefficients(tau_0_1, tau, beta,
coeffs_by_hand_1)
x_pts_1 = x_drive(t_pts, tau_0_1, tau, beta, N_max, coeffs_by_hand_1)
# now only tau is varied
tau = 1.5
coeffs_by_hand_2 = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
tau_0_2 = 1.
beta = 0.2
A_n_array_2, delta_n_array_2 = driven_coefficients(tau_0_2, tau, beta,
coeffs_by_hand_2)
x_pts_2 = x_drive(t_pts, tau_0_2, tau, beta, N_max, coeffs_by_hand_2)
# now only tau is varied
tau = 2.0
coeffs_by_hand_3 = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
tau_0_3 = 1.
beta = 0.2
A_n_array_3, delta_n_array_3 = driven_coefficients(tau_0_3, tau, beta,
coeffs_by_hand_3)
x_pts_3 = x_drive(t_pts, tau_0_3, tau, beta, N_max, coeffs_by_hand_3)
# now only tau is varied
tau = 2.5
coeffs_by_hand_4 = rectangular_coeffs_by_hand(N_max, tau, Delta_tau, f_max)
tau_0_4 = 1.
beta = 0.2
A_n_array_4, delta_n_array_4 = driven_coefficients(tau_0_4, tau, beta,
coeffs_by_hand_4)
x_pts_4 = x_drive(t_pts, tau_0_4, tau, beta, N_max, coeffs_by_hand_4)
How would you code the last cell differently so that you could just give an array of \(\tau\) values?
Compare coefficients for the cases from parts 1 to 4: large \(A_n\) coefficients when that component is in resonance.
print(' n case 1 case 2 case 3 case 4')
for n, (A_n_1, A_n_2, A_n_3, A_n_4) in enumerate(zip(A_n_array_1, A_n_array_2,
A_n_array_3, A_n_array_4)):
print(f'{n:2d} {A_n_1: .8f} {A_n_2: .8f} {A_n_3: .8f} {A_n_4: .8f}')
n case 1 case 2 case 3 case 4
0 0.00633257 0.00422172 0.00316629 0.00253303
1 0.17911224 0.01447102 0.00822068 0.00592958
2 0.00268521 0.00892472 0.08955612 0.01303477
3 0.00047497 0.00179014 0.00396133 0.00973767
4 0.00000000 0.00057109 0.00134261 0.00245255
5 -0.00009501 0.00015945 0.00056729 0.00107408
6 -0.00007678 0.00000000 0.00023749 0.00053672
7 -0.00003393 -0.00005543 0.00007835 0.00027238
8 -0.00000000 -0.00006360 0.00000000 0.00012820
9 0.00001584 -0.00005119 -0.00003562 0.00004629
10 0.00001629 -0.00003214 -0.00004751 0.00000000
11 0.00000864 -0.00001389 -0.00004630 -0.00002467
12 0.00000000 -0.00000000 -0.00003839 -0.00003583
13 -0.00000522 0.00000837 -0.00002778 -0.00003854
14 -0.00000591 0.00001158 -0.00001697 -0.00003608
15 -0.00000339 0.00001086 -0.00000745 -0.00003071
16 -0.00000000 0.00000774 -0.00000000 -0.00002399
17 0.00000233 0.00000372 0.00000509 -0.00001696
18 0.00000277 0.00000000 0.00000792 -0.00001036
19 0.00000167 -0.00000266 0.00000879 -0.00000462
Now for the graphs in Figure 5.25:
fig_2 = plt.figure(figsize=(10,5))
ax_1 = fig_2.add_subplot(2,2,1)
ax_1.plot(t_pts, x_pts_1, color='blue')
ax_1.set_xlim(0., 7.)
ax_1.set_ylim(-0.2, 0.2)
ax_1.set_xlabel('t')
ax_1.set_ylabel('x(t)')
ax_1.set_title(r'case (1): $\tau = 1$')
ax_2 = fig_2.add_subplot(2,2,2)
ax_2.plot(t_pts, x_pts_2, color='blue')
ax_2.set_xlim(0., 7.)
ax_2.set_ylim(-0.2, 0.2)
ax_2.set_xlabel('t')
ax_2.set_ylabel('x(t)')
ax_2.set_title(r'case (2): $\tau = 1.5$')
ax_3 = fig_2.add_subplot(2,2,3)
ax_3.plot(t_pts, x_pts_3, color='blue')
ax_3.set_xlim(0., 7.)
ax_3.set_ylim(-0.2, 0.2)
ax_3.set_xlabel('t')
ax_3.set_ylabel('x(t)')
ax_3.set_title(r'case (3): $\tau = 2.0$')
ax_4 = fig_2.add_subplot(2,2,4)
ax_4.plot(t_pts, x_pts_4, color='blue')
ax_4.set_xlim(0., 7.)
ax_4.set_ylim(-0.2, 0.2)
ax_4.set_xlabel('t')
ax_4.set_ylabel('x(t)')
ax_4.set_title(r'case (4): $\tau = 2.5$')
fig_2.tight_layout()
fig_2.savefig('Example_.png')
To add additional graphs with more \(\tau\) values, we should really define a function that produces the plot given a \(\tau\).