Fourier Series

(I’m continuing to use a Jupyter notebook, as in the last post, Least Squares.)

Sage has some rudimentary support for Fourier series, as part of the “piecewise-defined function” class, but it seems to be very slow and very flaky. Let’s implement our own. To make things run reasonably efficiently, we’re going to have Sage do numerical, rather than symbolic, integrals.

As a warm up, to integrate x^2 from 0 to 1, you type:
numerical_integral(x^2,0,1)
The first entry is the answer, while the second is an error estimate.

With this in hand, here’s a function to return a Fourier approximation:

def numerical_fourier(f, n):
    ans = 0
    for i in range(1,n+1):
        ans+= (numerical_integral(f(x)*sin(i*x), -pi,pi)[0])*sin(i*x)/float(pi)
        ans+= (numerical_integral(f(x)*cos(i*x), -pi,pi)[0])*cos(i*x)/float(pi)
    ans+= numerical_integral(f(x),-pi, pi)[0]/(2*float(pi))
    return ans

The inputs to this function are the function you want to approximate, and the number n of sine terms to use.

Let’s approximate the function f(x)=e^sin(x), which is a random-ish periodic function of period 2pi, by a Fourier series with 7 terms:
numerical_fourier(lambda x:exp(sin(x)), 3)


This “lambda x: exp(sin(x))” is a way of quickly defining and passing a function, called lambda calculus.

The answer is not so enlightening. It’s a bit better if we plot it:
plot(exp(sin(x)),(x,-2*pi,2*pi), color='red')
plot(numerical_fourier(lambda x:exp(sin(x)), 3), (x,-2*pi,2*pi))


Let’s plot several Fourier approximations on the same plot, to watch them converging:
plot(exp(sin(x)),(x,-2*pi,2*pi), color="red")+plot(numerical_fourier(lambda x:exp(sin(x)), 1), (x,-2*pi,2*pi), color='green')+plot(numerical_fourier(lambda x:exp(sin(x)), 2), (x,-2*pi,2*pi), color='blue')+plot(numerical_fourier(lambda x:exp(sin(x)), 3), (x,-2*pi,2*pi), color='orange')

It converges pretty quickly — it’s hard to even see the difference between orange (approximation) and red (exact).