#Script to illustrate how to compute a Fourier cosine series for a function f(x)
#defined on an interval 0 <= x <= L.

#Define function f(x)
var('f x');
f = function('f')(x)
f(x) = x - x^2/3;

#Interval length L and number of terms in Fourier cosine expansion
L = 2;
n = 5;

#Set up array to hold Fourier coefficients
a = vector(SR, [0..n]);

#Loop k = 0 to k = n, compute Fourier cosine coefficients of f(x)
for k in range(0,n+1):
    a[k] =  2/L*integral(f(x)*cos(k*pi*x/L),x,0,L);

#Form the Fourier cosine approximation
var('j');
fcos = a[0]/2 + sum(a[j]*cos(j*pi*x/L) for j in (1..n));

#Display result
print(fcos)

#Plot the function f(x) and Fourier cosine approximation
plt1 = plot(f(x),(x,0,L),color='blue');
plt2 = plot(fcos,(x,0,L),color='red');
show(plt1+plt2)