AMAT 415 — Matlab activity 2: Plotting Fourier series

In this activity, we will plot some Fourier series and see how well they represent the functions that give rise to them. More precisely, if f(t) is a function, we will compare its graph to the graph of its Fourier series

\displaystyle{\sum_{n=-\infty}^\infty \hat{f}(n)e^{int},\quad\text{where}\quad \hat{f}(t)=\frac{1}{2\pi}\int_{-\pi}^\pi f(t)e^{-int}dt}.

Let N be an integer. Enter N=10 to start with, say, and define NN=-N:N;. Since we can’t add up infinitely many terms of a series in practice, we’ll actually evaluate

\displaystyle{\sum_{n=-N}^N \hat{f}(n)e^{int}}.

The domain of our function f(t) is the interval [-\pi,\pi]. We’ll represent our domain as M+1 points spaced evenly in [-\pi,\pi]. Enter M=100 to start with and define the vector T=-pi:2*pi/M:pi;.

For each entry j between 1 and length(NN) and each k between 1 and length(T), store the value exp(i*NN(j)*T(k)) as the (j,k)-th entry of a matrix E. Letting E=exp(i*NN’*T); should do the trick.

The function f(t) still has not played a role. It enters now. We will start with a rectangular pulse function: f=T>0;. You can plot it — plot(x,f).

To compute the fourier coefficients we need to integrate. To keep things simple, let’s just use the trapezoidal rule for doing the numerical integration. That way we know what’s happening “under the hood.” See the Wikipedia article if you’ve forgotten how approximating integrals with the trapezoidal rule works. If you’d like to use more sophisticated algorithms, check out the command quad (for quadrature). Start by initializing a vector to store the Fourier coefficients in: coeffs=zeros(size(NN));. You can compute them using a for loop:

for m = 1:length(NN)
coeffs(m)=trapz(T,1/(2*pi)*f./E(m,:));
end

(I’ll buy a donut for the first person to show me how to compute coeffs in one line without using a loop.)

And now let’s sum the series: s = sum(diag(coeffs)*E);. (Why did that command sum the series?) Take a look at the entries of s. What do you notice about their imaginary parts? Now plot s and f together: plot(T,[f;s]);. Admire the picture.

Now play around:

  • Increase N and see the approximations improve.  Plot the approximations for multiple N on the same graph.
  • Use different functions f. Try f=T. Try some randomish polynomials of small degree.
  • Play with the parameter M.
  • For f=T.^2 or another polynomial with only even powers of T appearing,  max(abs(f-s)) versus N for N in some range.  (You’ll want to use a for loop here.)  What observations do you make about the rate of convergence?

I may ask you to hand in some of your pictures/code as part of Assignment 2, so save/print your two nicest pictures and the code you used to generate them.

Here’s a picture for f=T.^2-2*T.^3+4*T.^6-7*T.^7+0.5 with N=20 and M=1000.

trigpoly

This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *