import numpy as np
import matplotlib.pyplot as plt
nxs = 1000
x = np.linspace(-np.pi, np.pi, nxs)
dx = 2 * np.pi / (nxs - 1)
nterms = [2, 5, 10]
u = np.zeros((len(nterms), nxs))
ntermsmax = nterms[-1]
phi = np.zeros((ntermsmax + 1, nxs))
phi[0, :] = 1 / np.sqrt(2 * np.pi)
for m in range(1, ntermsmax):
phi[m, :] = 1 / np.sqrt(np.pi) * np.cos(m * x)
flist = np.exp(-8 * (x ** 2) / np.pi**2)
for i in range(len(nterms)):
for k in range(nterms[i] + 1):
ck = np.dot(flist, phi[k, :]) * dx
u[i, :] += ck * phi[k, :]
plt.plot(x, flist, '-', x, u[0, :], '--', x, u[1, :], ':', x, u[2, :], '-.')
plt.legend(['exp(-8x^2/pi^2)', 'K=2', 'K=5', 'K=10'])
plt.show(block=False)
table = np.column_stack((x, flist, u.transpose()))
np.savetxt('FourierGaussian.dat', table, fmt='%f')