Parallel Spectral Numerical Methods/Examples in Matlab
Contents
Examples in Matlab and Python[edit]
We now want to find approximate numerical solutions using Fourier spectral methods. In this section we focus primarily on the heat equation with periodic boundary conditions for . Many of the techniques used here will also work for more complicated partial differential equations for which separation of variables cannot be used directly.
1D Heat Equation[edit]
The 1D heat equation

(
is a well known second order PDE for which exact series solutions can be found using separation of variables. It arises in several contexts such as in predicting the temperature in a thin uniform cross section rod. The equation and its derivation can be found in introductory books on partial differential equations and calculus, for example , and , The constant is the thermal diffusivity and is temperature. We have already described how to solve the heat equation using separation of variables. Let us first discretize such that where . are uniformly spaced in . Let’s now take the FFT of both sides of the 1D heat equation to obtain
We then rewrite the spatial derivative using the equation
The subscript denotes the coefficient of the Fourier mode.
so that the partial differential equation now becomes a collection of independent ODEs. While we can solve these ODEs in time exactly, we will use techniques that will also allow us to obtain approximate solutions to PDEs we cannot solve exactly. We will discuss two methods for solving these ODEs, forward Euler and backward Euler.
Forward Euler[edit]
Using the forward Euler method in time, we obtain
All that is left is to take the IFFT of the computed solution after all timesteps are taken to transfer it back to real space. This is a linear PDE, so only one IFFT is needed at the end. We will later see that this is different for a nonlinear PDE. A Matlab implementation of this is in listing A .

( 
%Solving Heat Equation using pseudospectral and Forward Euler
%u_t= \alpha*u_xx
%BC= u(0)=0, u(2*pi)=0
%IC=sin(x)
clear all; clc;
%Grid
N = 64; %Number of steps
h = 2*pi/N; %step size
x = h*(1:N); %discretize xdirection
alpha = .5; %Thermal Diffusivity constant
t = 0;
dt = .001;
%Initial conditions
v = sin(x);
k=(1i*[0:N/21 0 N/2+1:1]);
k2=k.^2;
%Setting up Plot
tmax = 5; tplot = .1;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;
for i = 1:nplots
v_hat = fft(v); %Fourier Space
for n = 1:plotgap
v_hat = v_hat+dt*alpha*k2.*v_hat; %FE timestepping
end
v = real(ifft(v_hat)); %Back to real space
data(i+1,:) = v;
t=t+plotgap*dt;
tdata = [tdata; t]; %Time vector
end
%Plot using mesh
mesh(x,tdata,data), grid on,
view(60,55), xlabel x, ylabel t, zlabel u, zlabel u

( 
#!/usr/bin/env python
"""
Solving Heat Equation using pseudospectral 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 xdirection
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)
# Plot using mesh
xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('t')
plt.show()
Backward Euler[edit]
To derive this method, we start by applying the FFT and then perform timestepping using backward Euler. We then rewrite the implicit form into a form that gives the next iterate,
Figure i shows a numerical solution to the heat equation (Methods to obtain the exact solution can be found in, among other places, Boyce and DiPrima^{[1]}.) where obtained using the Matlab program in listing B .

( 
%Solving Heat Equation using pseudospectral methods with Backwards Euler:
%u_t= \alpha*u_xx
%BC = u(0)=0 and u(2*pi)=0 (Periodic)
%IC=sin(x)
clear all; clc;
%Grid
N = 64; h = 2*pi/N; x = h*(1:N);
% Initial conditions
v = sin(x);
alpha = .5;
t = 0;
dt = .001; %Timestep size
%(ik)^2 Vector
k=(1i*[0:N/21 0 N/2+1:1]);
k2=k.^2;
%Setting up Plot
tmax = 5; tplot = .1;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;
for i = 1:nplots
v_hat = fft(v); %Converts to fourier space
for n = 1:plotgap
v_hat = v_hat./(1dt*alpha*k2); %Backwards Euler timestepping
end
v = ifft(v_hat); %Converts back to real Space
data(i+1,:) = real(v); %Records data
t=t+plotgap*dt; %Records time
tdata = [tdata; t];
end
%Plot using mesh
mesh(x,tdata,data), grid on, %axis([1 2*pi 0 tmax 1 1]),
view(60,55), xlabel x, ylabel t, zlabel u, zlabel u,

( 
#!/usr/bin/env python
"""
Solving Heat Equation using pseudospectral methods with Backwards Euler:
u_t= \alpha*u_xx
BC = u(0)=0 and u(2*pi)=0 (Periodic)
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; h = 2*math.pi/N; x = [h*i for i in xrange(1,N+1)]
# Initial conditions
v = [math.sin(y) for y in x]
alpha = 0.5
t = 0
dt = .001 #Timestep size
# (ik)^2 Vector
I = complex(0,1)
k = numpy.array([I*n for n in range(0,N/2) + [0] + range(N/2+1,0)])
k2=k**2;
# Setting up Plot
tmax = 5.0; tplot = 0.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) # convert to fourier space
for n in xrange(plotgap):
v_hat = v_hat / (1dt*alpha*k2) # backward Euler timestepping
v = numpy.fft.ifft(v_hat) # convert back to real space
data[i+1,:] = numpy.real(v) # records data
t = t+plotgap*dt # records real time
tdata.append(t)
# Plot using mesh
xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('t')
plt.show()
Exercises[edit]

Write a program to solve the heat equation using the CrankNicolson method.

Solve the advection equation for with the initial data


for and for
up to a time . You can do this either by using separation of variables or by assuming that the solution is of the form and deducing what is in order to satisfy the initial conditions. In both cases please use the forward Euler, backward Euler and CrankNicolson timestepping schemes. After calculating the exact solution in each of these cases, examine how the maximum error at the final time depends on the timestep for each of these three methods.

Nonlinear Equations[edit]
The 1D AllenCahn Equation[edit]
So far we have dealt only with linear equations. Now it’s time for a nonlinear PDE. The AllenCahn equation models the separation of phases in a material. It was introduced by S. Allen and J. W. Cahn^{[2]} and is

(
where is a small but positive constant. The way to numerically solve this is similar to the method used for the heat equation, but there are some notable differences. The biggest difference is that FFT(), so the must be computed before taking the FFT. The FFT is a linear operation but cubing is nonlinear operation, so the order matters

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \frac{\partial \hat{u}_k}{\partial t} = \epsilon\frac{\partial^2 \hat{u}_k}{\partial x^2}+\hat{u}_k\widehat{u^3}_k. }
(
Next rewrite the first term on the right hand side, just like we did in the heat equation
In order to solve this numerically we are going to use a combination of implicit (backward Euler) and explicit (forward Euler) methods. We are going to skip forward Euler because it is unstable.
ImplicitExplicit Method[edit]
You might have already noticed that backward Euler is not going to work for the AllenCahn in its present state because of the nonlinear term. If you go to implement backward Euler you can see that you can’t factor out all of the . Luckily there is a simple intuitive way around this that isn’t detrimental to the accuracy of the solution. Write all the terms implicitly (backwards Euler) except for the nonlinear term which is expressed explicitly. Applying this to AllenCahn we find that ^{[3]}
Now we have a form that we can work with. We can set the initial conditions to be Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle u(x,0)=\frac{1}{4} \sin(x)} and plot the computed spacetime evolution calculated by the Matlab code in listing C . The computed result is in Fig. ii .

( 
%Solving 1D AllenCahn Eq using pseudospectral and Implicit/Explicit method
%u_t=u_{xx} + u  u^3
%where uu^3 is treated explicitly and u_{xx} is treated implicitly
%BC = u(0)=0, u(2*pi)=0 (Periodic)
%IC=.25*sin(x);
clear all; clc;
%Grid and Initial Data
N = 8000; h = 2*pi/N; x = h*(1:N); t = 0;
dt = .001; %timestep size
epsilon= .001;
%initial conditions
v = .25*sin(x);
%(ik) and (ik)^2 vectors
k=(1i*[0:N/21 0 N/2+1:1]);
k2=k.^2;
%setting up plot
tmax = 5; tplot = .2;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;
for i = 1:nplots
for n = 1:plotgap
v_hat = fft(v); %converts to Fourier space
vv = v.^3; %computes nonlinear term in real space
vv = fft(vv); %converts nonlinear term to Fourier space
v_hat = (v_hat*(1/dt+1)  vv)./(1/dtk2*epsilon); %Implicit/Explicit
v = ifft(v_hat); %Solution back to real space
end
data(i+1,:) = real(v); %Records data each "plotgap"
t=t+plotgap*dt; %Real time
tdata = [tdata; t];
end
mesh(x,tdata,data), grid on, axis([1 2*pi 0 tmax 1 1]),
view(60,55), xlabel x, ylabel t, zlabel u
The 2D AllenCahn Equation[edit]
Now we will look at the 2D form of the AllenCahn Equation, where satisfies

(
The convert it into Fourier space by taking the FFT of both sides

(
where and is to remind us that we take the FFT in respected directions. We will also define

(
The way to deal with the first two terms on the right hand side is to take the FFT in the xdirection and then take it in the ydirection. The order in which the FFT is done, first or first is not important. Some software libraries offer a two dimensional FFT. It usually depends on the equation being solved whether it is more efficient to use a multidimensional FFT or many one dimensional FFTs. Typically, it is easier to write a program which uses a multidimensional FFT, but in some situations this is not very efficient since one can immediately reuse data that has just been Fourier transformed.
ImplicitExplicit Method[edit]
In this method, the nonlinear term given in eq. 6 is calculated explicitly, while the rest of the terms will be written implicitly such that
we can then substitute in for

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \hat{u}_k^{n+1} = \frac{\hat{u}_k^n\left(\frac{1}{h}+1\right)\widehat{(u^n)^3}_k}{ \left(\epsilon(ik_x)^2\epsilon(ik_y)^2+\frac{1}{h}\right)}. }
(
The Matlab code used to generate Fig. iii is in listing D .

( 
%Solving 2D AllenCahn Eq using pseudospectral with Implicit/Explicit
%u_t= epsilon(u_{xx}+u_{yy}) + u  u^3
%where uu^3 is treated explicitly and epsilon(u_{xx} + u_{yy}) is treated implicitly
%BC = Periodic
%IC=v=sin(2*pi*x)+0.001*cos(16*pi*x;
clear all; clc;
%Grid
N = 256; h = 1/N; x = h*(1:N);
dt = .01;
%x and y meshgrid
y=x';
[xx,yy]=meshgrid(x,y);
%initial conditions
v=sin(2*pi*xx)+0.001*cos(16*pi*xx);
epsilon=.01;
%(ik) and (ik)^2 vectors in x and y direction
kx=(2*pi*1i*[0:N/21 0 N/2+1:1]);
ky=(2*pi*1i*[0:N/21 0 N/2+1:1]');
k2x=kx.^2;
k2y=ky.^2;
[kxx,kyy]=meshgrid(k2x,k2y);
for n = 1:500
v_nl=v.^3; %calculates nonlinear term in real space
%FFT for linear and nonlinear term
v_nl = fft2(v_nl);
v_hat=fft2(v);
vnew=(v_hat*(1+1/dt)v_nl)./ ...
((kxx+kyy)*epsilon+1/dt); %Implicit/Explicit timestepping
%converts to real space in xdirection
v=ifft2(vnew);
%Plots each timestep
mesh(v); title(['Time ',num2str(n)]); axis([0 N 0 N 1 1]);
xlabel x; ylabel y; zlabel u;
view(43,22); drawnow;
end

( 
#!/usr/bin/env python
"""
Solving 2D AllenCahn Eq using pseudospectral with Implicit/Explicit
u_t= epsilon(u_{xx}+u_{yy}) + u  u^3
where uu^3 is treated explicitly and u_{xx} and u_{yy} is treated implicitly
BC = Periodic
IC=v=sin(2*pi*x)+0.5*cos(4*pi*y)
"""
import math
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import time
plt.ion()
# Setup the grid
N = 64; h = 1.0/N;
x = [h*i for i in xrange(1,N+1)]
y = [h*i for i in xrange(1,N+1)]
dt = 0.05
xx,yy = (numpy.mat(A) for A in (numpy.meshgrid(x,y)))
# Initial Conditions
u = numpy.array(numpy.sin(2*math.pi*xx) + 0.5*numpy.cos(4*math.pi*yy), dtype=float)
epsilon = 0.01
# (ik) and (ik)^2 vectors in x and y direction
I = complex(0,1)
k_x = 2*numpy.pi*numpy.array([I*n for n in range(0,N/2) + [0] + range(N/2+1,0)])
k_y = k_x
kxx = numpy.zeros((N,N), dtype=complex)
kyy = numpy.zeros((N,N), dtype=complex)
for i in xrange(N):
for j in xrange(N):
kxx[i,j] = k_x[i]**2
kyy[i,j] = k_y[j]**2
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xx, yy, u,rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
v_hat = numpy.zeros((N,N), dtype=complex)
v_hat = numpy.fft.fft2(u)
for n in xrange(100):
# calculate nonlinear term in real space
v_nl = numpy.array(u**3, dtype=complex)
# FFT for nonlinear and linear terms
v_nl = numpy.fft.fft2(v_nl)
v_hat = (v_hat*(1+1/dt)  v_nl)
v_hat=v_hat/(1/dt  (kxx+kyy)*epsilon) # Implicit/Explicit timestepping
u = numpy.real(numpy.fft.ifft2(v_hat))
# Remove old plot before drawing
ax.collections.remove(surf)
surf = ax.plot_surface(xx, yy, u,rstride=1, cstride=1, cmap=cm.jet,
linewidth=0, antialiased=False)
plt.draw()
plt.show()
Exercises[edit]
Many of these exercises are taken from Uecker^{[4]}. Another introductory source of information on these equations is Trefethen and Embree^{[5]}.
 Burgers equation is given by: where and has periodic boundary conditions. Solve this equation using an implicitexplicit method. If you take to be small, ensure that a sufficient number of grid points are used to get the correct numerical solution. A simple way to check this is to keep increasing the number of grid points and checking that there is no change in the solution. Another way to check this is to calculate the Fourier coefficients and check that the highest ones decay to machine precision.
 The KuramotoSivashinsky equation is given by: where has periodic boundary conditions.
 What does this equation model and what type of behavior do you expect its solutions to have?
 Find numerical solutions to this equation using an implicitexplicit method.
 The 1D GrayScott equations are given by: where , , and are constants.
 What does this equation model and what type of behavior do you expect its solutions to have?
 Find numerical solutions to this equation using an implicitexplicit method. Try several different values of , , and and compare the resulting patterns to what you can find in the literature.
 The 2D SwiftHohenberg equation is given by:
 What does this equation model and what type of behavior do you expect its solutions to have?
 Find numerical solutions to this equation using an implicitexplicit method for several values of .
 The 2D GrayScott equations are given by: where , , and are constants.
 What does this equation model and what type of behavior do you expect its solutions to have?
 Find numerical solutions to this equation using an implicitexplicit method.
 The 2D Complex GinzburgLandau equation is given by: An introductory tutorial to this equation can be found at http://codeinthehole.com/static/tutorial/index.html
 What does this equation model and what type of behavior do you expect its solutions to have?
 Find numerical solutions to this equation using an implicitexplicit method for several values of and .
Notes[edit]
 ↑ Boyce and DiPrima (2010)
 ↑ Allen and Cahn (1979)
 ↑ Notice that when programming you are going to have to update the nonlinear term () each time you want to calculate the next timestep . The reason this is worth mentioning is because for each timestep you are going to have to go from real space to Fourier space to real space, then repeat. For, the heat equation you can perform any number of timesteps in Fourier space and only convert back when you record your data.
 ↑ Uecker (2009)
 ↑ Trefethen and Embree
References[edit]
Allen, S.M.; Cahn, J.W. (1979). "A microscopic theory for antiphase boundary motion and its applications to antiphase domain coarsening". Acta Metallurgica 27: 10851095.
DiPrima, R.C. (2010). Elementary Differential Equations and Boundary Value Problems. Wiley.
Trefethen, L.N., ed.,; Embree, K., ed.,, The (Unfninished) PDE coffee table book, http://people.maths.ox.ac.uk/trefethen/pdectb.html
Uecker, H.. A short ad hoc introduction to spectral methods for parabolic PDE and the NavierStokes equations. http://www.staff.unioldenburg.de/hannes.uecker/hfwebe.html.