\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts}

\newcommand{\x}{\boldsymbol{x}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}

\begin{document}

\title{Numerical Methods I -- Problem Set 8}
\author{Fall Semester 2005}
\date{Due November 14}

\maketitle

\noindent
In \emph{trigonometric interpolation} we use a truncated Fourier
series to interpolate a periodic function $u$ on the interval
$[0,2\pi]$.  Let $N$ be an even number of equidistant nodes $x_j=jh$,
where $h=2\pi/N$ and $j=0,\dots,N-1$.

The trigonometric interpolant of $u$ is the function
\begin{equation}
  v(x) = \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x} \,.
  \label{e.reconstruction}
\end{equation}
To determine the coefficients $c_k$, we impose the interpolation
condition $v(x_j)=u(x_j)$ for $j=0, \dots, N-1$, i.e.
\begin{equation}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j} = u(x_j) \,.
\end{equation}
We multiply this equation by $\e^{-\i m x_j}$ and sum over all $j$, so
that
\begin{equation}
  \sum_{j=0}^{N-1} \e^{-\i m x_j}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j} =
  \sum_{j=0}^{N-1} \e^{-\i m x_j} \, u(x_j) \,.
\end{equation}
Notice that
\begin{align}
  \sum_{j=0}^{N-1} \e^{-\i m x_j}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j}
  & = \sum_{k=-N/2}^{N/2-1} c_k \sum_{j=0}^{N-1} \e^{\i (k-m) x_j} \notag \\
  & = \sum_{k=-N/2}^{N/2-1} c_k \sum_{j=0}^{N-1} \e^{\i (k-m) h j} \notag \\
  & = \sum_{k=-N/2}^{N/2-1} c_k \cdot
  \begin{cases}
    N & \text{if } k=m \\
    \dfrac{1 - (\e^{\i (k-m) 2\pi/N})^N}{1 - \e^{\i (k-m) 2\pi/N}} = 0
    & \text{if } k \neq m 
  \end{cases} \notag \\
  & = N \, c_m \,.
\end{align}
We therefore obtain that
\begin{equation}
  c_m = \frac1N \sum_{j=0}^{N-1} \e^{-\i m x_j} \, u(x_j) \,.
  \label{e.dft}
\end{equation}

\paragraph{Remark:} The normalization convention is arbitrary.  Most
library routines, including those in \textsf{Octave}, have the factor
$1/N$ in the inverse transform \eqref{e.reconstruction} rather than in
the forward transform \eqref{e.dft}.  The convention above has the
advantage that the \emph{discrete Fourier transform} \eqref{e.dft} can
be seen as the Riemann sum approximation of the continuous Fourier
transform.  If both the forward and the inverse transform get a factor
of $1/\sqrt N$, the transform is unitary, which has certain
theoretical advantages, but is usually avoided in actual code.

\begin{enumerate}

\item \textbf{Project:} Write an \textsf{Octave} program which
computes the coefficients $c_m$ for $m=-\tfrac{N}2,\dots,\tfrac{N}2-1$
for a given function $u$, e.g.
\[
   u(x) = x \, (x-2 \, \pi) \, \e^{-x} \,.
\]
This operation can be done very efficiently via the \emph{Fast Fourier
Transform}, available in \textsl{Octave} as the function \texttt{fft}.
Compare your result with that produced by \texttt{fft} to find out how
the coefficients $c_m$ are laid out in \textsf{Octave}'s data
structure. 

\item \textbf{Project:} Write an \textsf{Octave} function
\texttt{tpolyval(c,x)}, which, in analogy with the built-in function
\texttt{polyval}, evaluates the trigonometric interpolant $v$ via
equation \eqref{e.reconstruction}.  Use the coefficient layout as
employed by \texttt{fft}.

\begin{enumerate}
\item Plot $u$ and $v$ for the given example.

\item What happens when $N=10$ and you interpolate $u(x)=\sin(6x)$?

\end{enumerate}
\end{enumerate}

\end{document}
