Dealing with boundary conditions using Fourier spectral methods

Question

I am currently working on a project where I need to use Fourier spectral methods to solve the KS equation. I found this code which is using the Fourier spectral methods to solve the classic 1D heat equation. Here's the code:

"""
Solving Heat Equation using pseudo-spectral and Forward Euler
u_t= \alpha*u_xx
BC= u(0)=0, u(2*pi)=0
IC=sin(x)
"""

import math
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator

# Grid
N = 64 # Number of steps
h = 2*math.pi/N # step size
x = h*numpy.arange(0,N) # discretize x-direction

alpha = 0.5 # Thermal Diffusivity constant
t = 0
dt = .001

# Initial conditions
v = numpy.sin(x)
I = complex(0,1)
k = numpy.array([I*y for y in range(0,N/2) + [0] + range(-N/2+1,0)])
k2=k**2;

# Setting up Plot
tmax = 5; tplot = .1;
plotgap = int(round(tplot/dt))
nplots = int(round(tmax/tplot))

data = numpy.zeros((nplots+1,N))
data[0,:] = v
tdata = [t]

for i in xrange(nplots):
    v_hat = numpy.fft.fft(v)

    for n in xrange(plotgap):
        v_hat = v_hat+dt*alpha*k2*v_hat # FE timestepping

    v = numpy.real(numpy.fft.ifft(v_hat)) # back to real space
    data[i+1,:] = v

    # real time vector
    t = t+plotgap*dt
    tdata.append(t)

I don't really understand how the code deals with the boundary conditions. I understand that the initial condition satisfies the boundary conditions. But in the code there is nothing about the boundary conditions.

How should I proceed to change the boundary conditions in this code? For exemple how would I compute periodic boundary conditions for the heat equation?

Answer 1

In this code, you are trying to solve this 1D heat transfer equation:

$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^{2} u}{\partial x^{2}}$$

If you use discrete Fourier transform (DFT) on $u(x,t)$ and discretized $x$ values of $x_{i}$ where $0 \leq i \leq N$ and $x \in [0,L]$ and $x_{i} = \frac{i}{N}L$ and $k_{i} = \frac{2 \pi i}{L}$, you have:

$$\tilde{u}(k_{m},t) = \sum_{n=0}^{N} u(x_{n},t) e^{-i x_{n} k_{m}}$$

and

$$u(x_{n},t) = \sum_{m = 0}^{N} \tilde{u}(k_{m},t) e^{i x_{n} k_{m}}$$

The first assumption here is that $u$ is periodic. Because:

$$u(x_{n}+L,t) = \sum_{m=0}^{N} \tilde{u}(k_{m},t) e^{i x_{n} k_{m}} e^{i L k_{m}} = \sum_{m=0}^{N} \tilde{u}(k_{m},t) e^{i x_{n} k_{m}} e^{i 2 \pi m} = u(x_{n},t)$$

So, you have this transformed governing equation as:

$$\frac{d \tilde{u} (k,t)}{dt} = -\alpha k^{2} \tilde{u}(k,t)$$

You can discretized this equation by using Forward Euler method as:

$$\tilde{u}(k,t+\Delta t) = \tilde{u}(k,t) + \Delta t (-\alpha k^{2} \tilde{u}(k,t))$$

You know that: $u(x,0) = \sin(x)$, so you have $\tilde{u}(k,0)$. So keep in mind in this scheme you don't see the boundary conditions explicitly cause they are hidden in the assumtpion of DFT that $u$ must be periodic.