# Import symbolics, plotting, linear algebra, Legendre
# polynomials, and numerical quadrature functions.
from pylab import plot, savefig, grid, legend, clf, pcolor, spy
from numpy import arange, zeros, sqrt
from math import (pi, factorial, sin, cos, log, exp, tan, atan, asin,
acos, sinh, cosh, tanh, asinh, acosh)
from scipy.linalg import solve
from scipy.special.orthogonal import legendre
from scipy.special.orthogonal import p_roots
from scipy.integrate import quadrature
# Define Fourier basis.
def fn(i, x):
if i == 0:
return 1/sqrt(2*pi)
if i % 2 == 0:
return sin((i // 2)*x)/sqrt(pi)
else:
return cos((i // 2 + 1)*x)/sqrt(pi)
# Evaluate at 'x' linear combination of fn(i)
# with the coefficients vector 'vec'
def lin_comb(x, vec):
val = 0
n = len(vec)
for i in range(0, n):
val = val + vec[i]*fn(i, x)
return val
# Integrate product of f(x) with fn(i)
def gauss_quad(f, fn, i, quad_x, quad_w):
val = 0
for j in range(0, len(quad_x)):
val += f(quad_x[j])*fn(i, quad_x[j])*quad_w[j]
return val
# Calculate L2 error norm
def calculate_L2_error(f, rhs, quad_x, quad_w):
err = 0
for j in range(0, len(quad_x)):
val = f(quad_x[j]) - lin_comb(quad_x[j], rhs)
err += val**2 * quad_w[j]
return sqrt(err)
# Calculate projection coefficients. Use the fact that functions
# 1/sqrt(2*pi), cos(x)/sqrt(pi), sin(x)/sqrt(pi), cos(2x)/sqrt(pi),
# sin(2x)/sqrt(pi) ... are orthogonal on (-pi,pi) in the L2-product
# (and thus no matrix inversion is necessary).
def calc_fourier_coeffs(f, n, num_gauss_points, print_rhs):
rhs = zeros(n) # zero vector
# fill the rhs vector
print "Preparing Gauss quadrature..."
roots, weights = p_roots(num_gauss_points)
print "Done."
print "num_gauss_points =", num_gauss_points
print "WARNING: You should experiment with num_gauss_points"
print " and increase it until results converge."
quad_x = [pi*xx for xx in roots]
quad_w = [pi*w for w in weights]
print "Calculating integrals..."
# calculating integrals of f times fn(i) in (-pi, pi)
for i in range(0, n):
rhs[i] = gauss_quad(f, fn, i, quad_x, quad_w)
print "Done."
# showing right-hand side rhs
if print_rhs: print "Fourier coefficients:", rhs
# calculate L2 error numerically
print "Calculating L2 error..."
err_L2 = calculate_L2_error(f, rhs, quad_x, quad_w)
print "L2 error:", err_L2
return rhs
# Plot the projection result, which is a linear combination of
# fn(i) with the coefficients from the vector vec.
def plot_result(vec, f, m):
print "Plotting result..."
h_ref = 2*pi/m # 2*pi is the length of the reference interval (-pi, pi)
x_array = arange(-pi, pi+h_ref, h_ref)
y_array = [lin_comb(xx, vec) for xx in x_array]
f_array = [f(xx) for xx in x_array]
e_array = [f(xx) - lin_comb(xx, vec) for xx in x_array]
clf()
label = "Fourier series F(x)"
plot(x_array, y_array, "b-", label=label)
label = "Function f(x)"
plot(x_array, f_array, "g-", label=label)
legend()
grid(True)
savefig("a.png")
clf()
label = "Error E(x) = f(x) - F(x)"
plot(x_array, e_array, "b-", label=label)
legend()
grid(True)
savefig("err.png")
print "Done."
# approximate error maximum (over plotting points)
print "Calculating max error..."
err_abs = [abs(f(xx) - lin_comb(xx, vec)) for xx in x_array]
err_max = max(err_abs)
print "Max error:", err_max
return
# Plot the basis.
def plot_basis(n, m):
print "Plotting basis..."
h_ref = 2*pi/m # 2*pi is the length of the reference interval (-pi, pi)
x_array = arange(-pi, pi+h_ref, h_ref)
for i in range(0, n):
clf()
label = "Fourier basis function fn%d(x)" % (i)
y_array = [fn(i, xx) for xx in x_array]
plot(x_array, y_array, "b-", label=label)
legend()
grid(True)
savefig("four%d.png" % (i))
print "Done."
return
# Construct and plot the Fourier approximation F(x) of length 'n'
# of a given function f(x) defined in interval (a,b).
def fourier(f, n, plot_basis_fn=False, print_rhs=True, m=500):
n = int(n)
# WARNING: Since the Fourier functions are not polynomials and nor
# is the function f(x) in general, quadrature of very high order
# may be necessary. The following number of Gauss points is chosen
# heuristically. You should experiment with it and increase it if needed.
num_gauss_points = 4*n + 20 # remember: quadrature order = 2*this -1
# Perform orthogonal L2-projection of f(x) onto the linear
# space generated by the functions fn(0), ..., fn(n-1)
vec = calc_fourier_coeffs(f, n, num_gauss_points, print_rhs)
# plot the approximation F(x) along with f(x) in interval (-pi, pi)
plot_result(vec, f, m)
# Plot functions fn(0), ..., fn(n-1) in interval (-pi, pi)
if plot_basis_fn: plot_basis(n, m)
return
