Parallel Spectral Numerical Methods/Examples in Matlab and Python

Examples in Matlab and Python

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 $x\in [0,2\pi )$. 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

The 1D heat equation

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

( 1)

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 $\alpha$ is the thermal diffusivity and $u(x,t)$ is temperature. We have already described how to solve the heat equation using separation of variables. Let us first discretize $x$ such that $x_{j}$ where $j=0,1,2,...,n$. $x_{j}$ are uniformly spaced in $[0,2\pi )$. Let’s now take the FFT of both sides of the 1D heat equation to obtain

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

We then rewrite the spatial derivative using the equation ${\text{FFT}}({\frac {\partial ^{\nu }u_{j}}{\partial x^{\nu }}})\equiv (ik)^{\nu }{\hat {u}}_{k}.$

The subscript $k$ denotes the coefficient of the $k^{th}$ Fourier mode.

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\alpha (ik)^{2}{\hat {u}}_{k},$

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

Using the forward Euler method in time, we obtain

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\alpha (ik)^{2}{\hat {u}}_{k}^{n}$

${\hat {u}}_{k}^{n+1}={\hat {u}}_{k}^{n}+\alpha h(ik)^{2}{\hat {u}}_{k}^{n}$

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 .

( A)

A Matlab program to solve the heat equation using forward Euler timestepping Code Download

%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)
clear all; clc;

%Grid
N = 64; %Number of steps
h = 2*pi/N; %step size
x = h*(1:N); %discretize x-direction

alpha = .5; %Thermal Diffusivity constant
t = 0;
dt = .001;

%Initial conditions
v = sin(x);
k=(1i*[0:N/2-1 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

( Ap)

A Python program to solve the heat equation using forward Euler time-stepping. Code download

#!/usr/bin/env python
"""
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)

# 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

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,

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\alpha (ik)^{2}{\hat {u}}_{k}$

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\alpha (ik)^{2}{\hat {u}}_{k}^{n+1}$

${\hat {u}}_{k}^{n+1}(1-\alpha h(ik)^{2})={\hat {u}}_{k}^{n}$

${\hat {u}}_{k}^{n+1}={\frac {{\hat {u}}_{k}^{n}}{(1-\alpha h(ik)^{2})}}.$ 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.) where $n=64$ obtained using the Matlab program in listing B .

( i)

A numerical solution to the heat equation, eq. 1 computed using the backward Euler method.

( B)

A Matlab program to solve the heat equation using backward Euler timestepping Code Download

%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/2-1 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./(1-dt*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,

( Bp)

A Python program to solve the heat equation using backward Euler time-stepping. Code download

#!/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 / (1-dt*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

1. Write a program to solve the heat equation using the Crank-Nicolson method.

2. Solve the advection equation $u_{t}=u_{x}$ for $x\in [0,2\pi )$ with the initial data

   1. $u(t=0,x)=\cos(x)$

   2. $u(t=0,x)=0$ for $x<\pi$ and $u(t=0,x)=1$ for $x\geq \pi$

   up to a time $T=1$. You can do this either by using separation of variables or by assuming that the solution is of the form $u(x,t)=f(x+t)$ and deducing what $f$ is in order to satisfy the initial conditions. In both cases please use the forward Euler, backward Euler and Crank-Nicolson 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

The 1D Allen-Cahn Equation

So far we have dealt only with linear equations. Now it’s time for a nonlinear PDE. The Allen-Cahn equation models the separation of phases in a material. It was introduced by S. Allen and J. W. Cahn and is

${\frac {\partial u}{\partial t}}=\epsilon {\frac {\partial ^{2}u}{\partial x^{2}}}+u-u^{3},$

( 2)

where $\epsilon$ 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($u^{3}$)$\neq ({\text{FFT}}(u))^{3}$, so the $u^{3}$ must be computed before taking the FFT. The FFT is a linear operation but cubing is non-linear operation, so the order matters

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon {\frac {\partial ^{2}{\hat {u}}_{k}}{\partial x^{2}}}+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}.$

( 3)

Next rewrite the first term on the right hand side, just like we did in the heat equation

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon (ik)^{2}{\hat {u}}_{k}+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}.$

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.

Implicit-Explicit Method

You might have already noticed that backward Euler is not going to work for the Allen-Cahn 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 ${\hat {u}}_{k}^{n+1}$. 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 Allen-Cahn we find that

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\epsilon (ik)^{2}{\hat {u}}_{k}^{n+1}+{\hat {u}}_{k}^{n}-{\widehat {(u^{n})^{3}}}_{k}$

${\hat {u}}_{k}^{n+1}\left(-\epsilon (ik)^{2}+{\frac {1}{h}}\right)={\frac {1}{h}}{\hat {u}}_{k}^{n}+{\hat {u}}_{k}^{n}-{\widehat {(u^{n})^{3}}}_{k}$

${\hat {u}}_{k}^{n+1}={\frac {{\hat {u}}_{k}^{n}({\frac {1}{h}}+1)-{\widehat {(u^{n})^{3}}}_{k}}{\left(-\epsilon (ik)^{2}+{\frac {1}{h}}\right)}}.$

Now we have a form that we can work with. We can set the initial conditions to be $u(x,0)={\frac {1}{4}}\sin(x)$ and plot the computed space-time evolution calculated by the Matlab code in listing C . The computed result is in Fig. ii .

( C)

A Matlab program to solve the 1D Allen-Cahn equation using implicit explicit timestepping Code download

%Solving 1D Allen-Cahn Eq using pseudo-spectral and Implicit/Explicit method
%u_t=u_{xx} + u - u^3
%where u-u^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/2-1 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/dt-k2*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

( ii)

A numerical solution to the 1D Allen-Cahn equation, eq. 2 , with $\epsilon =0.001$ and $u(x,t=0)=0.25\sin(x)$ computed using an implicit explicit method.

The 2D Allen-Cahn Equation

Now we will look at the 2D form of the Allen-Cahn Equation, where $u(x,y,t)$ satisfies

${\frac {\partial u}{\partial t}}=\epsilon \left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)+u-u^{3}.$

( 4)

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

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon \left({\frac {\partial ^{2}{\hat {u}}_{k}}{\partial x^{2}}}+{\frac {\partial ^{2}{\hat {u}}_{k}}{\partial y^{2}}}\right)+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}$

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon \left((ik_{x})^{2}{\hat {u}}_{k}+(ik_{y})^{2}{\hat {u}}_{k}\right)+{\hat {u}}_{k}-{\widehat {(u^{3})}}_{k}$

( 5)

where $k_{x}$ and $k_{y}$ is to remind us that we take the FFT in respected directions. We will also define

$f(u)\equiv u-u^{3}$

( 6)

The way to deal with the first two terms on the right hand side is to take the FFT in the x-direction and then take it in the y-direction. The order in which the FFT is done, $x$ first or $y$ 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.

Implicit-Explicit Method

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

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\epsilon \left((ik_{x})^{2}{\hat {u}}_{k}^{n+1}+(ik_{y})^{2}{\hat {u}}_{k}^{n+1}\right)+{\widehat {f(u^{n})}}_{k}$

${\hat {u}}_{k}^{n+1}\left(-\epsilon (ik_{x})^{2}-\epsilon (ik_{y})^{2}+{\frac {1}{h}}\right)={\frac {{\hat {u}}_{k}^{n}}{h}}+{\widehat {f(u^{n})}}_{k}$

${\hat {u}}_{k}^{n+1}={\frac {{\frac {{\hat {u}}_{k}^{n}}{h}}+{\widehat {f(u^{n})}}_{k}}{\left(-\epsilon (ik_{x})^{2}-\epsilon (ik_{y})^{2}+{\frac {1}{h}}\right)}}$

we can then substitute in for $f(u)$

${\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)}}.$

( 7)

( iii)

A numerical solution to the 2D Allen-Cahn equation, eq. 4 at time $t=500$ with $\epsilon =0.1$ and $u(x,y,t=0)=\sin(2\pi x)+0.001\cos(16\pi x)$ computed using an implicit explicit method.

The Matlab code used to generate Fig. iii is in listing D .

( D)

A Matlab program to solve the 2D Allen-Cahn equation using implicit explicit timestepping Code download

%Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicit
%u_t= epsilon(u_{xx}+u_{yy}) + u - u^3
%where u-u^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/2-1 0 -N/2+1:-1]);
ky=(2*pi*1i*[0:N/2-1 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 x-direction
    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

( Dp)

A Python program to solve the 2D Allen Cahn equation using implicit explicit time-stepping. Code download

#!/usr/bin/env python
"""
Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicit
u_t= epsilon(u_{xx}+u_{yy}) + u - u^3
where u-u^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

Many of these exercises are taken from Uecker. Another introductory source of information on these equations is Trefethen and Embree.

1. Burgers equation is given by: ${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial x^{2}}}-u{\frac {\partial u}{\partial x}}$ where $\nu \in \mathbb {R} ^{+}$ and $u$ has periodic boundary conditions. Solve this equation using an implicit-explicit method. If you take $\nu$ 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.
2. The Kuramoto-Sivashinsky equation is given by: ${\frac {\partial u}{\partial t}}=-{\frac {\partial ^{2}u}{\partial x^{2}}}-{\frac {\partial ^{4}u}{\partial x^{4}}}-u{\frac {\partial u}{\partial x}}$ where $u$ 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 implicit-explicit method.
3. The 1D Gray-Scott equations are given by: ${\frac {\partial u}{\partial t}}=d_{1}{\frac {\partial ^{2}u}{\partial x^{2}}}-uv^{2}+f(1-u),$ ${\frac {\partial v}{\partial t}}=d_{2}{\frac {\partial ^{2}v}{\partial x^{2}}}+uv^{2}-(f+k)v$ where $d_{1}$, $d_{2}$, $f$ and $k$ 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 implicit-explicit method. Try several different values of $d_{1}$, $d_{2}$, $f$ and $k$ and compare the resulting patterns to what you can find in the literature.
4. The 2D Swift-Hohenberg equation is given by: ${\frac {\partial u}{\partial t}}=-\Delta ^{2}u+2\Delta u+(\alpha -1)u-u^{3},$
   • 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 implicit-explicit method for several values of $\alpha$.
5. The 2D Gray-Scott equations are given by: ${\frac {\partial u}{\partial t}}=d_{1}\Delta u-uv^{2}+f(1-u),$ ${\frac {\partial v}{\partial t}}=d_{2}\Delta v+uv^{2}-(f+k)v$ where $d_{1}$, $d_{2}$, $f$ and $k$ 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 implicit-explicit method.
6. The 2D Complex Ginzburg-Landau equation is given by: ${\frac {\partial A}{\partial t}}=A+(1+i\alpha )\Delta A-(1+i\beta )|A|^{2}A.$ 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 implicit-explicit method for several values of $\alpha$ and $\beta$.

References

Allen, S.M.; Cahn, J.W. (1979). "A microscopic theory for antiphase boundary motion and its applications to antiphase domain coarsening". Acta Metallurgica. 27: 1085–1095.

Boyce, W.E.; 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{{citation}}: CS1 maint: extra punctuation (link)

Uecker, H. A short ad hoc introduction to spectral methods for parabolic PDE and the Navier-Stokes equations.