Guide to writing a C++ program for computing the Fourier series

How to write a C++ program for computing the Fourier series

In this blog, we are diving into computational physics. And learn how to write a C++ program for the Fourier series 

Here's a step-by-step guide to writing a C++ program for computing the Fourier series of a periodic function. We'll assume the function is defined for one period [−π,π] for simplicity.

---

Step 1: Understand the Fourier Series

The Fourier series is represented as: $f(x)=a_0+{\sum }_{n=1}^{\mathrm{\infty }}(a_n\cos (nx)+b_n\sin (nx)) where:$

• $a_0=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }f(x)dx$
• 
  
• $a_n=\frac{1}{\pi }{\int }_{-\pi }^{\pi }f(x)\cos (nx)dx$
• 
  
• $b_n=\frac{1}{\pi }{\int }_{-\pi }^{\pi }f(x)\sin (nx)dx$
  

---

Step 2: Plan the Program Structure

1. Define the function $f(x)$
2. Use numerical integration (like the trapezoidal rule) to compute the coefficients $a_0,a_n,b_{\mathrm{n.}}$
3. Compute the Fourier series up to a finite number of terms $N$.
4. Output or visualize the result.

Step 3: Write the Program

Here's the code:

cpp
#include <iostream>
#include <cmath>
#include <vector>

// Function to integrate (e.g., f(x) = x^2)
double f(double x) {
    return x * x;
}

// Numerical integration using the trapezoidal rule
double integrate(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n; // Step size
    double sum = 0.5 * (func(a) + func(b)); // Boundary values

    for (int i = 1; i < n; i++) {
        sum += func(a + i * h);
    }

    return sum * h;
}

// Compute Fourier coefficients a0, an, and bn
void computeCoefficients(double &a0, std::vector<double> &an, std::vector<double> &bn, int N, int numSteps) {
    a0 = integrate(f, -M_PI, M_PI, numSteps) / (2 * M_PI);

    for (int n = 1; n <= N; n++) {
        auto cosTerm = [n](double x) { return f(x) * cos(n * x); };
        auto sinTerm = [n](double x) { return f(x) * sin(n * x); };

        an[n - 1] = integrate(cosTerm, -M_PI, M_PI, numSteps) / M_PI;
        bn[n - 1] = integrate(sinTerm, -M_PI, M_PI, numSteps) / M_PI;
    }
}

// Approximate the Fourier series
double fourierSeries(double x, double a0, const std::vector<double> &an, const std::vector<double> &bn, int N) {
    double result = a0;

    for (int n = 1; n <= N; n++) {
        result += an[n - 1] * cos(n * x) + bn[n - 1] * sin(n * x);
    }

    return result;
}

int main() {
    int N = 10;         // Number of Fourier terms
    int numSteps = 1000; // Steps for numerical integration

    double a0;
    std::vector<double> an(N), bn(N);

    // Compute coefficients
    computeCoefficients(a0, an, bn, N, numSteps);

    // Test the Fourier series
    for (double x = -M_PI; x <= M_PI; x += 0.1) {
        double approx = fourierSeries(x, a0, an, bn, N);
        std::cout << "x = " << x << ", f(x) = " << f(x) << ", Fourier Approximation = " << approx << std::endl;
    }

    return 0;
}

---

Step 4: Explanation

1. f(x) Function: Represents the function you want to approximate (e.g., $x^2$).
2. Integration: The integrate function uses the trapezoidal rule to compute integrals.
3. Coefficients: The computeCoefficients function calculates $a_0$, $a_n$, and $b_n$.
4. Fourier Approximation: The fourierSeries function reconstructs the Fourier series up to $N$ terms.
5. Main Function: Calculates and prints the approximation for $f(x)$ at various points.

---

Step 5: Compile and Run

1. Save the code as fourier_series.cpp.
2. Compile it using g++ fourier_series.cpp -o fourier_series -lm.
3. Run the program with ./fourier_series.

---

Step 6: Visualize (Optional)

Use tools like GNUplot or libraries like Matplotlib (Python) to visualize the Fourier series approximation against the original function.

For more blogs like this stay connected to Physics Blog