Scientific Computing, Spring 1999

Assignment 7.

Given April 20, due end of finals week.

Objective: To code the FFT and see an application

This assignment makes use of a Fourier sine series, which is closely related to exponential Fourier series we have worked with. Because e^{it} = cos(t) + i sin(t), the exponential Fourier coefficients,

f̂_{α} = \frac{1}{2π} ∫_{0}^{2π} e^{−iαt}f(t)dt ,

can be written as the sum of a cosine coefficient and a sine coefficient:

f̂_{α} = a_{α} − ib_{α} ,

where

a_{α} = \frac{1}{2π} ∫_{0}^{2π} cos(αt)f(t)dt ,

b_{α} = \frac{1}{2π} ∫_{0}^{2π} sin(αt)f(t)dt .

In a similar way, we can define a discrete cosine transform and a discrete sin transform from the discrete (exponential) Fourier transform.

The sine and cosine transform can be useful in certain circumstances, particularly when the function, f, has special properties. In particular, the Fourier sine series is useful if f is real and skew symmetric about t = π. This means that f(π + t) = −f(π − t) for any t. Such functions may seem odd and unlikely to arise in practice. They do arise in the following way. Suppose that f(t) is initially defined for t only in the range from 0 to π, and that f is continuous on that range, and that f(0) = 0 and f(π) = 0. In solving differential equations, these would be called “Dirichlet boundary conditions”. A skew symmetric periodic function that is continuous satisfies these conditions automatically. Conversely, for any continuous function defined on the interval [0, π] with Dirichlet boundary values, we can define a periodic skew symmetric function, defined for all t, with the same values in the original [0, π] interval. First extend f to the complementary interval [π, 2π] by skew symmetry: f(π+t) = −f(π−t) for 0 ≤ t ≤ π. Then define f for t outside the fundamental period [0, 2π] by f(t+2πn) = f(t) for t ∈ [0, 2π] and n any integer. Of course, we never actually do this extension in the computer. Instead, we think of doing it to help us derive properties of the Fourier sine series using what we already know about the exponential Fourier series.

We can make a discrete version of this by thinking at first of the numbers f_{k} as being samples of a function: f_{k} = f(t_{k}). If we have N points t_{k} = kΔt, Δt = 2π/N, then the midpoint, π will be one of the t_{k} only if N is even. Therefore, suppose N = 2(M + 1) and that we have M points uniformly spaced inside the interval [0, π]:

t_{k} = kΔt , Δt = \frac{π}{M + 1} , for k = 1, …, M.  (1)

We can (mentally) extend the samples f_{k}, defined for k = 1, …, M to be defined for k = 0, …, N − 1 by skew symmetry. First define f_{0} = 0 (consistent with the Dirichlet boundary condition). Then define f_{k+M+1} = −f_{k} for k = 0, …, M.

The exponential Fourier coefficients of the extended f give the Fourier sine coefficients of the original f. In the continuous case,

f̂_{α} = −ib_{α} , where b_{α} = \frac{1}{π} ∫_{0}^{π} sin(αt)f(t)dt .  (2)

and

f(t) = 2 ∑_{α=1}^{∞} b_{α} sin(αt) .  (3)

The factor of 2 in (2) comes from the fact that the contributions from α and −α are equal, and we have kept only positive α terms. All the cosine coefficients, a_{α} are zero in this case, because the integrals

∫_{0}^{2π} cos(αt)f(t)dt

are zero.

1. Derive the discrete version of the Fourier sine series. This starts with M samples (real numbers) , f_{1}, …, f_{M}, and produces M coefficients b̃_{1}, …, b̃_{M}. You can find these formulae by using the DFT applied to the skew symmetric extended f.

2. Suppose that M = 2^{L} − 1. Derive a fast algorithm for computing the discrete sine coefficients that uses only real arithmetic. This will be called the FST (fast sine transform)

3. Assuming that the f_{k} are samples of a smooth, skew symmetric function, find an FST based algorithm that gives a spectrally accurate estimate of f′′(t_{k}) for k = 1, …, M. Call this operation A. That is, there is an M × M real matrix so that if Af = g, then g_{k} is the estimate of f′′(t_{k}) from the samples f_{j}, j = 1, …, M. Write a program to calculate Af using teh FST. Check numerically that the matrix A is symmetric. Check that it produces spectrally accurate results. For the latter, you will have to be careful in choosing f.

4. Write a program to solve the partial differential equation

\frac{∂u}{∂t} = \frac{1}{2} \frac{∂^{2}u}{∂x^{2}} + u^{2} ,  (4)

with boundary conditions

u(x = 0, t) = 0 , u(x = π, t) = 0 ,  (5)

and intitial conditions

u(x, t = 0) = z sin(x) .  (6)

Find an approximate solution by replacing u by a discrete approximation u_{k}(t) ≈ u(x_{k}, t), where x_{k} = kΔx, Δx = π/(M + 1), and k = 1, …, M. Approximate the equation (4) by

u̇_{k}(t) = \frac{1}{2} (Au(t))_{k} + u_{k}^{2} .  (7)

This is a system of M ordinary differential equations. Compute the operation of A using the algorithm from part 3. Solve these ordinary differential equations using the fourth order Runge Kutta method. You must take the time step, Δt, to be very small for large N. At a minimum, you must take Δt < \frac{1}{4}Δx^{2}, or the numerical solution will “blow up”^{1}. Do a convergence study to verify that your method has spectral accuracy. Compute solutions for various values of z. For small z the solutions tranquilly go to zero as t → ∞. For larger z, the solution develops a spike near x = π which leads to infinite^{2} u at x = π for some finite t. Observe how the computational method breaks down as the spike forms.

1This blow up has nothing to do with physical blowup of corectly computed solutions. It is a numerical artifact assiciated with errors in the discrete approximation.

2This blow up is real, and very hard to compute, as you will see.