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

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&rsquo;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 Matlab program to solve the heat equation using forward Euler timestepping   Code Download

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

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



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

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

Exercises
  Write a program to solve the heat equation using the Crank-Nicolson method.   Solve the advection equation $$u_t=u_x$$ for $$x\in[0,2\pi)$$ with the initial data   $$u(t=0,x)=\cos(x)$$   $$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. 

The 1D Allen-Cahn Equation
So far we have dealt only with linear equations. Now it&rsquo;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

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

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&rsquo;t factor out all of the $$\hat{u}_k^{n+1}$$. Luckily there is a simple intuitive way around this that isn&rsquo;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 $$. The computed result is in Fig. $$.

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



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

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

where $$k_x$$ and $$k_y$$ 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 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. $$ 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)$$



The Matlab code used to generate Fig. $$ is in listing $$.

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

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

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.
 * 3) * What does this equation model and what type of behavior do you expect its solutions to have?
 * 4) * Find numerical solutions to this equation using an implicit-explicit method.
 * 5) 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.
 * 6) * What does this equation model and what type of behavior do you expect its solutions to have?
 * 7) * 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.
 * 8) The 2D Swift-Hohenberg equation is given by: $$\frac{\partial u}{\partial t}=-\Delta^2u+2\Delta u+ (\alpha-1)u-u^3,$$
 * 9) * What does this equation model and what type of behavior do you expect its solutions to have?
 * 10) * Find numerical solutions to this equation using an implicit-explicit method for several values of $$\alpha$$.
 * 11) 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.
 * 12) * What does this equation model and what type of behavior do you expect its solutions to have?
 * 13) * Find numerical solutions to this equation using an implicit-explicit method.
 * 14) The 2D Complex Ginzburg-Landau equation is given by: $$\frac{\partial A}{\partial t}=A+(1+i\alpha)\Delta A- (1+i\beta)|A|^2A.$$ An introductory tutorial to this equation can be found at http://codeinthehole.com/static/tutorial/index.html
 * 15) * What does this equation model and what type of behavior do you expect its solutions to have?
 * 16) * Find numerical solutions to this equation using an implicit-explicit method for several values of $$\alpha$$ and $$\beta$$.