CODEWITHSUNDEEP
pythonnumerical-integrationcontinuous-fouriersimpsons-rule
Aran G
Aran G

Reputation: 135

Plotting Fourier Series coefficients in Python using Simpson's Rule

I want to 1. express Simpson's Rule as a general function for integration in python and 2. use it to compute and plot the Fourier Series coefficients of the function sin(t).

I've stolen and adapted this code for Simpson's Rule, which seems to work fine for integrating simple functions such as e^t, sin(t) or poly

Given period period, the Fourier Series coefficients are computed as:

a_k

b_k

where k = 1,2,3,...

I am having difficulty figuring out how to express ak. I'm aware that a_k=0 since this function is odd, but I would like to be able to compute it in general for other functions.

Here's my attempt so far:

import matplotlib.pyplot as plt
from numpy import *

def f(t):
    k = 1
    for k in range(1,10000): #to give some representation of k's span
        k += 1
    return sin(t)*sin(k*t)

def trapezoid(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for j in range(1, n):
        s += f(a + j*h)
    s += f(b)/2.0
    return s * h

print trapezoid(f, 0, 2*pi, 100)

This doesn't give the correct answer of 0 at all since it increases as k increases and I'm sure I'm approaching it with tunnel vision in terms of the for loop. My difficulty in particular is with stating the function so that k is read as k = 1,2,3,...

The problem I've been given unfortunately doesn't specify what the coefficients are to be plotted against, but I am assuming it's meant to be against k.

Upvotes: 0

Views: 1539

Answers (1)

Mateen Ulhaq
Mateen Ulhaq

Reputation: 27201

Here's one way to do it, if you want to run your own integration or fourier coefficient determination instead of using numpy or scipy's built in methods:

import numpy as np

def integrate(f, a, b, n):
    t = np.linspace(a, b, n)
    return (b - a) * np.sum(f(t)) / n

def a_k(f, k):
    def ker(t): return f(t) * np.cos(k * t)
    return integrate(ker, 0, 2*np.pi, 2**10+1) / np.pi

def b_k(f, k):
    def ker(t): return f(t) * np.sin(k * t)
    return integrate(ker, 0, 2*np.pi, 2**10+1) / np.pi

print(b_k(np.sin, 0))

This gives the result

0.0

On a side note, trapezoid integration is not very useful for uniform time intervals. But if you desire:

def trap_integrate(f, a, b, n):
    t = np.linspace(a, b, n)
    f_t = f(t)
    dt = t[1:] - t[:-1]
    f_ab = f_t[:-1] + f_t[1:]
    return 0.5 * np.sum(dt * f_ab)

There's also np.trapz if you want to use pre-builtin functionality. Similarly, there's also scipy.integrate.trapz

Upvotes: 1

Related Questions