#!/usr/bin/python # ================================================================ # John Kerl # kerl.john.r@gmail.com # 2006-03-03 # # This is an example of truncated Fourier series in one variable. # ================================================================ # This software is released under the terms of the GNU GPL. # Please see LICENSE.txt in the same directory as this file. # ================================================================ # ---------------------------------------------------------------- from __future__ import division # 1/2 = 0.5, not 0. from math import * import sys # ---------------------------------------------------------------- # Numerical quadrature for computing Fourier coefficients. def fourcoeff(func, k, xlo, xhi, nx, dx): sum = 0.0 x = xlo for i in range(0, nx): expnikx = cos(k*x) - 1j*sin(k*x) sum += expnikx * func(x) * dx x += dx return sum / 2.0 / pi # ---------------------------------------------------------------- # Function to be approximated via truncated Fourier series def f(x): #return x #return pi - abs(x) #return cos(x) #return x^2 #return pi**2 - x**2 #return cos(1.7*x) #return cos(1.7*x+0.4) #return x**3 #return cos(4*x) + cos(5*x) + cos(6*x) + 1.2 * cos(7*x) #return exp(x) #return 1 / (1 + x**2) #return cosh(x) #return sinh(x) #return exp(-pi*x**2) #a = pi/5 a = pi/2 if (abs(x) <= a): return 1.0 else: return 0.0 # ================================================================ # Start of main program xlo = -pi; xhi = pi; nx = 400; dx = (xhi - xlo) / (1.0 * nx); kmax = 20; # ---------------------------------------------------------------- # First, compute the Fourier coefficients. # Index the ck array from 0 to 2*kmax+1. ck = range(0, 2*kmax+1) for k in range (-kmax, kmax+1): ck[k+kmax] = fourcoeff(f, k, xlo, xhi, nx, dx) if 1: print "Fourier coefficients:" for k in range (-kmax, kmax+1): print "%3d: %11.7f %11.7f" % (k, ck[k+kmax].real, ck[k+kmax].imag) print # ---------------------------------------------------------------- # Next, tabulate the approximations. x = xlo for i in range (0, nx): print "%11.7f" % (x), print "%11.7f" % (f(x)), y = ck[kmax+0] print "%11.7f" % (y.real), for k in range (1, kmax): cposk = ck[kmax+k] expikx = cos(k*x) + 1j*sin(k*x) y += cposk * expikx cnegk = ck[kmax-k] expnikx = cos(k*x) - 1j*sin(k*x) y += cnegk * expnikx print "%11.7f" % (y.real), print x += dx