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:



