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

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)])

# 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

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?

  • 2
    $\begingroup$ BCs are hidden in the assumption of using FFT here. When you use FFT on original governing equation (here it is 1D heat transfer equation) and then just use the initial conditon and inverse FFT to get the solution back, you simply assumed that your domain (1D line here) is just periodic. $\endgroup$ – Alone Programmer Jan 25 '20 at 15:51
  • 1
    $\begingroup$ @alone programmer please make that an answer ... $\endgroup$ – Brian Borchers Jan 25 '20 at 22:25

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}}$$


$$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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy