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

The 1D heat equation

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

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

We then rewrite the spatial derivative using the equation${\displaystyle \text{FFT}(\frac{{\mathrm{\partial }}^{\nu }{u}_j}{\mathrm{\partial }{x}^{\nu }})\equiv (ik{)}^{\nu }{\hat{u}}_k.}$

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

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\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.

Using the forward Euler method in time, we obtain

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

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

%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)clearall;clc;%GridN=64;%Number of stepsh=2*pi/N;%step sizex=h*(1:N);%discretize x-directionalpha=.5;%Thermal Diffusivity constantt=0;dt=.001;%Initial conditionsv=sin(x);k=(1i*[0:N/2-10-N/2+1:-1]);k2=k.^2;%Setting up Plottmax=5;tplot=.1;plotgap=round(tplot/dt);nplots=round(tmax/tplot);data=[v;zeros(nplots,N)];tdata=t;fori=1:nplotsv_hat=fft(v);%Fourier Spaceforn=1:plotgapv_hat=v_hat+dt*alpha*k2.*v_hat;%FE timesteppingendv=real(ifft(v_hat));%Back to real spacedata(i+1,:)=v;t=t+plotgap*dt;tdata=[tdata;t];%Time vectorend%Plot using meshmesh(x,tdata,data),gridon,view(-60,55),xlabelx,ylabelt,zlabelu,zlabelu

#!/usr/bin/env python"""Solving Heat Equation using pseudo-spectral and Forward Euleru_t= \alpha*u_xxBC= u(0)=0, u(2*pi)=0IC=sin(x)"""importmathimportnumpyimportmatplotlib.pyplotaspltfrommpl_toolkits.mplot3dimportAxes3Dfrommatplotlibimportcmfrommatplotlib.tickerimportLinearLocator# GridN=64# Number of stepsh=2*math.pi/N# step sizex=h*numpy.arange(0,N)# discretize x-directionalpha=0.5# Thermal Diffusivity constantt=0dt=.001# Initial conditionsv=numpy.sin(x)I=complex(0,1)k=numpy.array([I*yforyinrange(0,N/2)+[0]+range(-N/2+1,0)])k2=k**2;# Setting up Plottmax=5;tplot=.1;plotgap=int(round(tplot/dt))nplots=int(round(tmax/tplot))data=numpy.zeros((nplots+1,N))data[0,:]=vtdata=[t]foriinxrange(nplots):v_hat=numpy.fft.fft(v)forninxrange(plotgap):v_hat=v_hat+dt*alpha*k2*v_hat# FE timesteppingv=numpy.real(numpy.fft.ifft(v_hat))# back to real spacedata[i+1,:]=v# real time vectort=t+plotgap*dttdata.append(t)# Plot using meshxx,tt=(numpy.mat(A)forAin(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()

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,

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\partial }t}=\alpha (ik{)}^2{\hat{u}}_k}$

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

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

${\displaystyle {\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^[1].) where${\displaystyle n=64}$ 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)clearall;clc;%GridN=64;h=2*pi/N;x=h*(1:N);% Initial conditionsv=sin(x);alpha=.5;t=0;dt=.001;%Timestep size%(ik)^2 Vectork=(1i*[0:N/2-10-N/2+1:-1]);k2=k.^2;%Setting up Plottmax=5;tplot=.1;plotgap=round(tplot/dt);nplots=round(tmax/tplot);data=[v;zeros(nplots,N)];tdata=t;fori=1:nplotsv_hat=fft(v);%Converts to fourier spaceforn=1:plotgapv_hat=v_hat./(1-dt*alpha*k2);%Backwards Euler timesteppingendv=ifft(v_hat);%Converts back to real Spacedata(i+1,:)=real(v);%Records datat=t+plotgap*dt;%Records timetdata=[tdata;t];end%Plot using meshmesh(x,tdata,data),gridon,%axis([-1 2*pi 0 tmax -1 1]),view(-60,55),xlabelx,ylabelt,zlabelu,zlabelu,

#!/usr/bin/env python"""Solving Heat Equation using pseudospectral methods with Backwards Euler:u_t= \alpha*u_xxBC = u(0)=0 and u(2*pi)=0 (Periodic)IC=sin(x)"""importmathimportnumpyimportmatplotlib.pyplotaspltfrommpl_toolkits.mplot3dimportAxes3Dfrommatplotlibimportcmfrommatplotlib.tickerimportLinearLocator# GridN=64;h=2*math.pi/N;x=[h*iforiinxrange(1,N+1)]# Initial conditionsv=[math.sin(y)foryinx]alpha=0.5t=0dt=.001#Timestep size# (ik)^2 VectorI=complex(0,1)k=numpy.array([I*nforninrange(0,N/2)+[0]+range(-N/2+1,0)])k2=k**2;# Setting up Plottmax=5.0;tplot=0.1plotgap=int(round(tplot/dt))nplots=int(round(tmax/tplot))data=numpy.zeros((nplots+1,N))data[0,:]=vtdata=[t]foriinxrange(nplots):v_hat=numpy.fft.fft(v)# convert to fourier spaceforninxrange(plotgap):v_hat=v_hat/(1-dt*alpha*k2)# backward Euler timesteppingv=numpy.fft.ifft(v_hat)# convert back to real spacedata[i+1,:]=numpy.real(v)# records datat=t+plotgap*dt# records real timetdata.append(t)# Plot using meshxx,tt=(numpy.mat(A)forAin(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()

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

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

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

   2. ${\displaystyle u(t=0,x)=0}$ for and${\displaystyle u(t=0,x)=1}$ for${\displaystyle x\ge \pi }$

   up to a time${\displaystyle T=1}$. You can do this either by using separation of variables or by assuming that the solution is of the form${\displaystyle u(x,t)=f(x+t)}$ and deducing what${\displaystyle 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.

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

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

( 2)

where${\displaystyle ϵ}$ 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(${\displaystyle u^3}$)${\displaystyle \ne (\text{FFT}(u){)}^3}$, so the${\displaystyle 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

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\partial }t}=ϵ\frac{{\mathrm{\partial }}^2{\hat{u}}_k}{\mathrm{\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

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\partial }t}=ϵ(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.

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${\displaystyle {\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^[3]

${\displaystyle \frac{{\hat{u}}_k^{n+1}-{\hat{u}}_k^n}{h}=ϵ(ik{)}^2{\hat{u}}_k^{n+1}+{\hat{u}}_k^n-{\widehat{(u^n{)}^3}}_k}$

${\displaystyle {\hat{u}}_k^{n+1}\left(-ϵ(ik{)}^2+\frac{1}{h}\right)=\frac{1}{h}{\hat{u}}_k^n+{\hat{u}}_k^n-{\widehat{(u^n{)}^3}}_k}$

${\displaystyle {\hat{u}}_k^{n+1}=\frac{{\hat{u}}_k^n(\frac{1}{h}+1)-{\widehat{(u^n{)}^3}}_k}{\left(-ϵ(ik{)}^2+\frac{1}{h}\right)}.}$

Now we have a form that we can work with. We can set the initial conditions to be${\displaystyle 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.

%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);clearall;clc;%Grid and Initial DataN=8000;h=2*pi/N;x=h*(1:N);t=0;dt=.001;%timestep sizeepsilon=.001;%initial conditionsv=.25*sin(x);%(ik) and (ik)^2 vectorsk=(1i*[0:N/2-10-N/2+1:-1]);k2=k.^2;%setting up plottmax=5;tplot=.2;plotgap=round(tplot/dt);nplots=round(tmax/tplot);data=[v;zeros(nplots,N)];tdata=t;fori=1:nplotsforn=1:plotgapv_hat=fft(v);%converts to Fourier spacevv=v.^3;%computes nonlinear term in real spacevv=fft(vv);%converts nonlinear term to Fourier spacev_hat=(v_hat*(1/dt+1)-vv)./(1/dt-k2*epsilon);%Implicit/Explicitv=ifft(v_hat);%Solution back to real spaceenddata(i+1,:)=real(v);%Records data each "plotgap"t=t+plotgap*dt;%Real timetdata=[tdata;t];endmesh(x,tdata,data),gridon,axis([-12*pi0tmax-11]),view(-60,55),xlabelx,ylabelt,zlabelu

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

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

( 4)

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

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\partial }t}=ϵ\left(\frac{{\mathrm{\partial }}^2{\hat{u}}_k}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{\hat{u}}_k}{\mathrm{\partial }{y}^2}\right)+{\hat{u}}_k-{\widehat{u^3}}_k}$

${\displaystyle \frac{\mathrm{\partial }{\hat{u}}_k}{\mathrm{\partial }t}=ϵ\left((i{k}_x{)}^2{\hat{u}}_k+(i{k}_y{)}^2{\hat{u}}_k\right)+{\hat{u}}_k-{\widehat{(u^3)}}_k}$

( 5)

where${\displaystyle k_x}$ and${\displaystyle k_y}$ is to remind us that we take the FFT in respected directions. We will also define

${\displaystyle 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,${\displaystyle x}$ first or${\displaystyle 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.

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

${\displaystyle \frac{{\hat{u}}_k^{n+1}-{\hat{u}}_k^n}{h}=ϵ\left((i{k}_x{)}^2{\hat{u}}_k^{n+1}+(i{k}_y{)}^2{\hat{u}}_k^{n+1}\right)+{\widehat{f(u^n)}}_k}$

${\displaystyle {\hat{u}}_k^{n+1}\left(-ϵ(i{k}_x{)}^2-ϵ(i{k}_y{)}^2+\frac{1}{h}\right)=\frac{{\hat{u}}_k^n}{h}+{\widehat{f(u^n)}}_k}$

${\displaystyle {\hat{u}}_k^{n+1}=\frac{\frac{{\hat{u}}_k^n}{h}+{\widehat{f(u^n)}}_k}{\left(-ϵ(i{k}_x{)}^2-ϵ(i{k}_y{)}^2+\frac{1}{h}\right)}}$

we can then substitute in for${\displaystyle f(u)}$

${\displaystyle {\hat{u}}_k^{n+1}=\frac{{\hat{u}}_k^n\left(\frac{1}{h}+1\right)-{\widehat{(u^n{)}^3}}_k}{\left(-ϵ(i{k}_x{)}^2-ϵ(i{k}_y{)}^2+\frac{1}{h}\right)}.}$

( 7)

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

%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;clearall;clc;%GridN=256;h=1/N;x=h*(1:N);dt=.01;%x and y meshgridy=x';[xx,yy]=meshgrid(x,y);%initial conditionsv=sin(2*pi*xx)+0.001*cos(16*pi*xx);epsilon=.01;%(ik) and (ik)^2 vectors in x and y directionkx=(2*pi*1i*[0:N/2-10-N/2+1:-1]);ky=(2*pi*1i*[0:N/2-10-N/2+1:-1]');k2x=kx.^2;k2y=ky.^2;[kxx,kyy]=meshgrid(k2x,k2y);forn=1:500v_nl=v.^3;%calculates nonlinear term in real space%FFT for linear and nonlinear termv_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-directionv=ifft2(vnew);%Plots each timestepmesh(v);title(['Time ',num2str(n)]);axis([0N0N-11]);xlabelx;ylabely;zlabelu;view(43,22);drawnow;end

#!/usr/bin/env python"""Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicitu_t= epsilon(u_{xx}+u_{yy}) + u - u^3where u-u^3 is treated explicitly and u_{xx} and u_{yy} is treated implicitlyBC = PeriodicIC=v=sin(2*pi*x)+0.5*cos(4*pi*y)"""importmathimportnumpyimportmatplotlib.pyplotaspltfrommpl_toolkits.mplot3dimportAxes3Dfrommatplotlibimportcmfrommatplotlib.tickerimportLinearLocatorimporttimeplt.ion()# Setup the gridN=64;h=1.0/N;x=[h*iforiinxrange(1,N+1)]y=[h*iforiinxrange(1,N+1)]dt=0.05xx,yy=(numpy.mat(A)forAin(numpy.meshgrid(x,y)))# Initial Conditionsu=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 directionI=complex(0,1)k_x=2*numpy.pi*numpy.array([I*nforninrange(0,N/2)+[0]+range(-N/2+1,0)])k_y=k_xkxx=numpy.zeros((N,N),dtype=complex)kyy=numpy.zeros((N,N),dtype=complex)foriinxrange(N):forjinxrange(N):kxx[i,j]=k_x[i]**2kyy[i,j]=k_y[j]**2fig=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)forninxrange(100):# calculate nonlinear term in real spacev_nl=numpy.array(u**3,dtype=complex)# FFT for nonlinear and linear termsv_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 timesteppingu=numpy.real(numpy.fft.ifft2(v_hat))# Remove old plot before drawingax.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()

Many of these exercises are taken from Uecker^[4]. Another introductory source of information on these equations is Trefethen and Embree^[5].

1. Burgers equation is given by:${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}=\nu \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}-u\frac{\mathrm{\partial }u}{\mathrm{\partial }x}}$ where${\displaystyle \nu \in {\mathbb{R}}^{+}}$ and${\displaystyle u}$ has periodic boundary conditions. Solve this equation using an implicit-explicit method. If you take${\displaystyle \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:${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}=-\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^4}-u\frac{\mathrm{\partial }u}{\mathrm{\partial }x}}$ where${\displaystyle 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:${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}=d_1\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}-u{v}^2+f(1-u),}$${\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}=d_2\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{x}^2}+u{v}^2-(f+k)v}$ where${\displaystyle d_1}$,${\displaystyle d_2}$,${\displaystyle f}$ and${\displaystyle 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${\displaystyle d_1}$,${\displaystyle d_2}$,${\displaystyle f}$ and${\displaystyle k}$ and compare the resulting patterns to what you can find in the literature.
4. The 2D Swift-Hohenberg equation is given by:${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}=-{\mathrm{\Delta }}^{2}u+2\mathrm{\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${\displaystyle \alpha }$.
5. The 2D Gray-Scott equations are given by:${\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}=d_1\mathrm{\Delta }u-u{v}^2+f(1-u),}$${\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }t}=d_2\mathrm{\Delta }v+u{v}^2-(f+k)v}$ where${\displaystyle d_1}$,${\displaystyle d_2}$,${\displaystyle f}$ and${\displaystyle 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:${\displaystyle \frac{\mathrm{\partial }A}{\mathrm{\partial }t}=A+(1+i\alpha )\mathrm{\Delta }A-(1+i\beta )|A{|}^{2}A.}$ An introductory tutorial to this equation can be found athttp://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${\displaystyle \alpha }$ and${\displaystyle \beta }$.

1. ↑Boyce and DiPrima (2010)
2. ↑Allen and Cahn (1979)
3. ↑Notice that when programming you are going to have to update the nonlinear term (${\displaystyle u^3}$) each time you want to calculate the next timestep${\displaystyle n+1}$. 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.
4. ↑Uecker (2009)
5. ↑Trefethen and Embree

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.

• Book:Parallel Spectral Numerical Methods

• CS1 maint: extra punctuation