Solving for coefficients of complex Fourier Series

Question

I am trying to prove the formula for the coefficients of the complex form Fourier Series. The context is that I want to be able solve for the constants using numerical integration, like Romberg's method. I have an old professor's proof to guide me, but I keep getting something different at one step.

Update 1: I figured out how the $m = n$ case plays into the proof. but still have one inconsistency at the end.

Update 2: Despite the lack of answers, I decided to try to implement this in python. The code is at the very bottom of the question.

\begin{align} f(x) =& \sum_{-\infty}^{\infty}c_n \exp(\frac{n \pi x}{L}i) \\\ f(x)\exp(-\frac{m \pi x}{L}i) =& \sum_{-\infty}^{\infty}c_n\exp(\frac{(n - m) \pi x}{L}i) \\\ \int_{-L}^{L} f(x)\exp(-\frac{m \pi x}{L}i) =& \sum_{-\infty}^{\infty} c_n\int_{-L}^{L} \exp(\frac{(n - m) \pi x}{L}i) \\\ =& \sum_{-\infty}^{\infty} c_n\int_{-L}^{L}\big[\cos(\frac{(n-m) \pi x}{L}) + i\sin(\frac{(n-m) \pi x }{L})\big]dx \\\ =& 2 \sum_{-\infty}^{\infty} c_n\int_{0}^{L}\cos(\frac{(n-m) \pi x}{L}) dx\\\ =& 2 \sum \frac{c_n L}{(n-m)\pi} \bigg[\sin(\frac{(n-m) \pi x}{L})\bigg]_0^L \\\ =&2 \sum \frac{c_n L}{(n-m)\pi} \bigg[\sin((n-m) \pi)\bigg] = 0 \\\ \text{if } m \not = n: \text{ } =& 2 \sum c_n \bigg[\sin(\frac{(n-m) \pi x}{L})\bigg]_0^L =0 \\\ \text{if } m = n: \text{ } =& \sum c_n \int_{-L}^{L}dx\\\ \implies \int_{-L}^{L} f(x)\exp(-\frac{m \pi x}{L}i) =& \sum c_n \delta_{n,m} \int_{-L}^{L}dx \\\ =& c_m \cdot 2L \\\ \implies c_m =& \frac{1}{2L}\int_{-L}^{L} f(x)\exp(-\frac{m \pi x}{L}i) \end{align}

So, if we want to compute $c_m$ with the composite trapezoid rule (for fixed step size):

\begin{align} \int_a^b g(x) dx \approx h\bigg[\frac{g(x_0) + g(x_N)}{2} + \sum_{j = 1}^{N - 1}g(x_j)\bigg] \end{align}

But here, we know that our fourier transform is periodic on $[-L, L]$ and therefore $f(x_0) = f(x_N)$, which gives us:

\begin{align} \int_a^b g(x) dx \approx h\bigg[\sum_{j = 0}^{N - 1}g(x_j)\bigg] \end{align}

And this is where I depart from the instructor's notes. My work follows:

Since $h = \frac{x_N - x_0}{N} =\frac{2L}{N}$ \begin{align} c_m =& \frac{h}{2L} \sum_{j = 0}^{N-1}f(x_j)\exp(\frac{im \pi x_j}{L}) \\\ =& \frac{1}{N}\sum_{j = 0}^{N-1}f(x_j)\exp(\frac{im \pi x_j}{L}) \end{align}

However, what my professor has: $$\sum_{j = 0}^{N-1}f(x_j)\exp(\frac{im \pi x_j}{L}) $$

What am I missing here? Is my form somehow equivalent, or did I fudge a step?

Here is the code implementation of the math above:

import matplotlib.pyplot as plt
import numpy as np


def coefficients(fn, dx, m, L):
    """
    Calculate the complex form fourier series coefficients for the first M
    waves.

    :param fn: function to sample
    :param dx: sampling frequency
    :param m: number of waves to compute
    :param L: We are solving on the interval [-L, L]
    :return: an array containing M Fourier coefficients c_m
    """

    N = 2*L / dx
    coeffs = np.zeros(m, dtype=np.complex_)
    xk = np.arange(-L, L + dx, dx)

    # Calculate the coefficients for each wave
    for mi in range(m):
        coeffs[mi] = 1/N * sum(fn(xk)*np.exp(-1j * mi * np.pi * xk / L))

    return coeffs


def fourier_graph(range, L, c_coef, function=None, plot=True, err_plot=False):
    """
    Given a range to plot and an array of complex fourier series coefficients,
    this function plots the representation.


    :param range: the x-axis values to plot
    :param c_coef: the complex fourier coefficients, calculated by coefficients()
    :param plot: Default True. Plot the fourier representation
    :param function: For calculating relative error, provide function definition
    :param err_plot: relative error plotted. requires a function to compare solution to
    :return: the fourier series values for the given range
    """
    # Number of coefficients to sum over
    w = len(c_coef)

    # Initialize solution array
    s = np.zeros(len(range))
    for i, ix in enumerate(range):
        for iw in np.arange(w):
            s[i] += c_coef[iw] * np.exp(1j * iw * np.pi * ix / L)

    # If a plot is desired:
    if plot:
        plt.suptitle("Fourier Series Plot")
        plt.xlabel(r"$t$")
        plt.ylabel(r"$f(x)$")
        plt.plot(range, s, label="Fourier Series")

        if err_plot:
            plt.plot(range, function(range), label="Actual Solution")
            plt.legend()

        plt.show()

    # If error plot is desired:
    if err_plot:
        err = abs(function(range) - s) / function(range)
        plt.suptitle("Plot of Relative Error")
        plt.xlabel("Steps")
        plt.ylabel("Relative Error")
        plt.plot(range, err)
        plt.show()

    return s


if __name__ == '__main__':

    # Assuming the interval [-l, l] apply discrete fourier transform:

    # number of waves to sum
    wvs = 50

    # step size for calculating c_m coefficients (trap rule)
    deltax = .025 * np.pi

    # length of interval for Fourier Series is 2*l
    l = 2 * np.pi

    c_m = coefficients(np.exp, deltax, wvs, l)

    # The x range we would like to interpolate function values
    x = np.arange(0, l, .01)
    sol = fourier_graph(x, l, c_m, np.exp, err_plot=True)

And here are the graphs that are generated:

Answer 1

Big thanks to @user753642 for spotting my mistakes over on the stackoverflow network:

I was computing the $c_n$ coefficients from 0 to $m$, where m is the number of wave functions in the sum. But by definitions the coefficients look like:

$$ c_m = \frac{1}{2L}\sum_{-\infty}^{\infty}c_n\delta_{n,m}\int_{-L}^Ldx = \frac{1}{2L}\int_{-L}^L f(x) \exp(\frac{-im\pi x}{L}) $$

So when I was reconstructing my function with the Fourier Series with $c_m$ with $n = 0 \dots m$, I should have been going from $n =-m/2 \dots m/2$.

This fixed everything, and allowed me to drop the weird N/2 factor that I was using.