Using the solutions of an ODE to construct a Fourier series

179 Views Asked by At

Consider the following ODE $$y''(t)+\omega_n^2y(t)=0,$$ with periodic boundary conditions $$y(t)=y(t+T).$$ This is an eigenvalue problem, where $\omega_n$ is the eigenvalue. Non-trivial solutions exist provided $$\omega_n=\frac{2\pi n}{T},$$ where $n$ is an integer. The corresponding eigenfunctions are $$y_n(t)=C_n\exp(i\omega_n t).$$ I would like to create a fourier series using solutions to the above equation such that $$g(t)=\sum_{n=-\infty}^\infty y_n(t).$$ I have tried to create a piece of python code that will do this for the case where $$g(t) = \exp(2i \pi t).$$ It seems to generate $g(t)$ okay but its prediction for $g'(t)$ is not good and I would like to know why. Here is the code:

# Program to produce a fourier series from the solutions of
# y'' + omega^2 y = 0
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import complex_ode

def f(t, y):
    # y[0] = y
    # y[1] = dy/dt
    # f[0] = dy/dt
    # f[1] = d^2y/dt^2
    return [y[1], -omega ** 2 * y[0]]

# Setup t grid
nt = 128
t_min = -0.5
t_max =  0.5
dt = (t_max - t_min) / (nt - 1)
t = np.linspace(t_min, t_max, nt)


# Calculate discrete fourier transform of g(t)
g   = np.exp(2j * np.pi * t)
g_t = 2j * np.pi * np.exp(2j * np.pi * t) # _t denotes d/dt
g_fft   = np.fft.fft(g)   / nt
g_t_fft = np.fft.fft(g_t) / nt


# Calculate required omega values
omega_array = 2 * np.pi * np.fft.fftfreq(nt, d = dt)

# Initialise to zero array
y = np.zeros((nt, 2), dtype = np.complex64)

# Calculate fourier series
i = -1
for omega in omega_array:
     i += 1

     y_dummy = np.zeros((nt, 2), dtype = np.complex64)

     # Initialise at t = t_min
     y0 = np.array([g_fft[i], g_t_fft[i]], dtype = np.complex64)

     # Setup ode
     r = complex_ode(f)
     r.set_initial_value(y0, t_min)

     # Initialise u_dummy
     y_dummy[0,:] = r.y

     # Iterate in t
     it = 0
     while r.successful() and r.t < t_max - dt / 2:
         it += 1
         y_dummy[it,:] = r.integrate(r.t+dt)

     y += y_dummy

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, np.real(y[:,0]))
ax.set_xlabel('t')
ax.set_title('Real($y$)')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, np.imag(y[:,0]))
ax.set_xlabel('t')
ax.set_title('Imag($y$)')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, np.real(y[:,1]))
ax.set_xlabel('t')
ax.set_title('Real($dy/dt$)')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, np.imag(y[:,1]))
ax.set_xlabel('t')
ax.set_title('Imag($dy/dt$)') 

Here are the outputted figures:

Real(y) Imag(y) Real(dy/dt) Imag(dy/dt)