# 5.1. Pseudo Spectral Method 算法说明¶

## 5.1.1. Introducing remarks¶

• the pseudo-spectral(PS) methods are methods to solve partial differential equantions(PDE)
• they originate roughly in 1970
• the PS methods have successfully been applied to
• fluid dynamics(turbulence modeling, weather predictions)
• non-linear waves
• seismic modeling
• MHD
• ...

## 5.1.2. Basic principles of the pseudo-spectral method¶

• the ‘pseudo-spectral’ in the method refers to the spatial part of a PDE
• example: a apatial PDE:
$\begin{split}L u(x)=s(x), x\in V \\ b.c.:f(u(y))=0, y\in \partial V\end{split}$

L: a spatial differential operator(e.g. $$L=\partial_{xx}+\partial_{yy}$$ , etc. )

• wanted: numerical solution $$u^N(x)$$ such that the residual R:
$R(x):=Lu^N(x)-s(x)$

is small-but how do we difine the smallness?

• general procedure:
1. choose a finite set of trial functions(expansion functions) $$\phi_j,j=0,\dots , N-1$$ and expand $$u^N$$ in these functions
$u^N(x)=\sum_{j=0}^{N-1}\hat{u_j}\phi_j(x)$
1. choose a set of test functions $$\chi_n,n=0,1,2,\dots, N-1$$ and demand that
$(\chi_n,R)=0 \ \ for\ n=0,1\dots N-1 (scalar product)$
• ‘spectral methods’ means that the trial functions $$\phi_n$$ form a basis for a certain space of global, smooth functions (e.g. Fourier polynomials)(global: extending over the whole spatial domain of interest)
• there are various spectral methods, classified according to the test functions $$\chi_n$$ :
Galerkin method, tau method, collocation or pseudo-spectral method
• collocation or pseudo-spectral method:
$\chi_n(x) = \delta(x-x_n),$

where the $$x_n (n=0,1,\dots,N-1)$$ are special points, the collocation points

• the smallness condition for the residual becomes:
$\sum_{j=0}^{N-1}\hat{u_j}L\phi_j(x_n)-s(x_n)=0, n = 0,,1,2,\dots,N-1$

N equations to determine the unknown N coeffients $$\hat{u_j}$$

• remark: the solution at the collocation points is exact, in between them we interpolate the solution
• what trial functions to choose?
1. periodic b.c.: trigonometric functions (Fourier series)
2. non-periodic b.c.: orthogonal polynomials(main candidate: Chebyshev polynomials)

## 5.1.3. Pseudo-spectral methods and Fourier transforms¶

### 5.1.3.1. Comparison to analytical Fourier method¶

• in
$\begin{split}Lu(x)=s(x),\ \ x\in V\\ b.c.:f(u(y))=0,\ y\in \partial V\end{split}$

assume 1-D, and e.g. $$L=partial_{zz}$$ , $$\partial_{zz}u=s(x)$$

• Fourier transform:
$(-ik_z)^2\hat{u}(k_z)=\hat{s}(k_z)$
• in principle, we want to do this numerically, but we have to make sure about a few points...

### 5.1.3.2. pseudo-spectral method, the Fourier case¶

• the aim is to find the expansion coefficients $$\hat{u}_j$$ such that the residual
$\sum_{j=0}^{N-1} \hat{u}_j L \phi_j(x_n)-s(x_n)=0,\ n=0,1,2,\dots,N-1$

or

$\sum_{j=0}^{N-1} \hat{u}_j L e^{-ik_jx_n}-s(x_n)=0,\ n=0,1,2,\dots,N-1$

vanishes. If L is linear, then $$Le^{-ik_jx_n}=h(k_j)e^{-ik_jx_n}$$

$\sum_{j=0}^{N-1} \hat{u}_j h(k_j)e^{-ik_jx_n}-s(x_n)=0,\ n=0,1,2,\dots,N-1$

the ‘trick’ is to choose (turning now to 1D for simplicity)

$x_n=n\Delta, n=0,1,2,\dots N-1$

and

$\begin{split}k_j= 2\pi j/(N \Delta), j=-N/2,\dots,N/2\\ \Delta:\textrm{spatial resotuion}\end{split}$
• $$z_n$$ and $$k_j$$ are equi-spaced, and the condition on the residual becomces:
$\sum_{j=0}^{N-1} \hat{u}_j h(k_j)e^{-i2\pi j n/N}-s(x_n)=0,\ n=0,1,2,\dots,N-1$
• we define the discrete Foureier Transform DFT as
$\hat{u}_j=\frac{1}{\sqrt{N}} \sum_{n=0}^{N-1}u_ne^{2 \pi i nj/N}$
• with $$u_n=u(x_n)$$ , and the inverse $$DFT^{-1}$$ as
$u_n=\frac{1}{\sqrt{N}} \sum_{n=-N/2}^{N/2}\hat{u}_je^{-2 \pi i nj/N}$
• it can be shown that with the specific choice of $$k_j$$ and $$x_n$$
$\frac{1}{\sqrt{N}} \sum_{n=0}^{N-1}e^{2\pi i j'/N}e^{-2 \pi i n j/N}=\delta(j-j')$

[algebraic proof, using $$\sum_{n=0}^{N-1}q^n=(1-q^N)/(1-q)$$ ]

so that

$u=DFT^{-1}(DFT(u))$

(but just and only at the collocation points, actually $$u_n=DFT^{-1}(DFT(u_n))$$ )

• then the condition on the residual
$\sum_{j=0}^{N-1} \hat{u}_j h(k_j)e^{-i2\pi j n/N}-s(x_n)=0,\ n=0,1,2,\dots,N-1$

can thus, using the DFT, be written as:

$DFT^{-1}[\hat{u}_j h(k_j)](x_n)-s(x_n)=0$

and, on applying DFT,

$\hat{u}_j h(k_j)(x_n)-DFT[s(x_n)](k_j)=0$

-> we can maipulate our equations numerically with the DFT analogously as we do treat equations analytically with the Fourier transform

Remark:
• $$x_n$$ and $$k_j$$ are equi-spaced only for trigonometric polynomials, every set of expansion functions has its own characteristic distribution of collocation points-equi-distribution is an exception (Chebychev, Legendre polynomials etc)
• the sets {$$u_n$$} and {$$u_j^*$$} are completely equivalent, they contain the same information

## 5.1.4. Time-dependent problems¶

• example: diffusion equation in 1D:
$\begin{split}\partial_tu=v\partial_{zz}u\\ u(z,0)=u_0(z)\end{split}$
• we consider the euqation only at the collocation points $$z_n=n\Delta,n=0,1,\dots,N-1$$ writing symbolically
$\begin{split}\partial_tu_n=v\partial_{zz}u_n\\ u(z_n,0)=u_0(z_n)\end{split}$

apply a spatial DFT:

$\begin{split}\partial_t\hat{u}_j=-v4\pi^2(j/N)^2\hat{u}_j\\ \hat{u}(k_j,0)=\hat{u_0}(k_j)\end{split}$

where $$j=-N/2,\dots,N/2$$

we have a set of N ODEs the temporal integration is done in Fourier space

### 5.1.4.1. Temporal integration¶

• The idea is to move the initial condition to Fourier space, and to do the temporal integration in Fourier space, since there we have ODEs
• since we ahve a set of ODEs, in principle every numerical scheme for integrating ODEs can be applied
• often good is Runge-Kutta 4th order, adaptive step-size
• 4th order RK:
$\begin{split}du/dt&=F(u,t) \ \ \ (u has N components)\\ u^{n+1}&=u^n+1/6(r_1+2r_2+2r_3+r_4)\\ r_1&=\Delta t F(u^n,t_n)\\ r_2&=\Delta t F(u^n+1/2r_1,t_n+1/2\Delta t)\\ r_3&=\Delta t F(u^n+1/2r_2,t_n+1/2\Delta t)\\ r_4&=\Delta t F(u^n+r_3,t_n+\Delta t)\\\end{split}$
• adaptive step-size:(for efficiency of the code)

advance $$\Delta t$$ , and also $$\Delta t/2+\Delta t/2$$,compare the results with prescribed accuracy, depending on the result make $$\Delta t$$ smaller or larger.

## 5.1.5. how to treat non-linearities¶

• assume there is a term $$\rho(z)u(z)$$ in the original PDE
• we are working in F-space, using DFT, so at a given time we have available $$\rho^*_j$$ and $$u^*_j$$
• $$\rho u$$ corresponds to a convolution in F-space, but convolutions are expensive(CPU time) and must be avoided (~:math:N^2 )
• the procedure to calculate $$(\rho u)^*_j$$ is as follow(~ $$N log_2N$$ ):
1. given at time t are $$\rho^*_j$$ and $$u^*_j$$
2. calculate $$\rho_n=DFT^{-1}(\rho^*_j)$$ and $$u_n=DFT^{-1}(u^*_j)$$
3. multiply and store $$w_n=\rho_nu_n$$
4. move $$w_n$$ to F-space, $$w^*_k=DFT(w_n)$$
5. use $$w^*_j\ for \ (\rho u)^*_j$$

## 5.1.6. Aliasing–de-aliansing¶

• the Foureier modes used are:
$\begin{split}e^{ik_jz_n}\\ \textrm{with wave-vectors} k_j=2\pi j/N \Delta, j=-N/2,\dots, N/2\\ \textrm{ and grid-points} z_n=n\Delta, n=0,1,\dots,N-1,\\ \textrm{i.e. the modes are } e^{2 \pi i j n/N}\end{split}$
• at the grid points $$x_n,e^{2\pi i n j /N}$$ equals:
$e^{2 \pi i(j+lN)n/Nf}, l=\dots,-2,-1,0,1,2,\dots$

this implies that modes with:

$k=2\pi (j+lN)/(N\Delta )$

contribute to the DFT as if they had

$k=2\pi j /N \Delta$

i.e. high k modes alias/bias the amplitude a lower k modes!

• example: for $$\Delta =1,N=8$$ , our wave vectors are
$k_j=-\pi,-3\pi/4,\dots,-\pi/4,0,\pi/4,\dots,3\pi/4,\pi$

now e.g. to $$k=\pi/4$$ also the modes $$k=9\pi/4,17\pi/4 \dots$$ etc. contribute! i.e. modes outside the k-range we model bias the modeled k-range

example: grid of N=8 points, $$\Delta =1$$ :

$$sin(z \pi/4)$$ and $$sin(z 9 \pi/4)$$ appear as being the same function when sampled.

First consequence of the aliasing effect:

prescribed functions such as initial conditions $$u(z,t=0)$$ or source functions $$s(z,t)$$ are best provided as superpositions fo the explicitly available modes,

$u(z_n,t=0)=\sum_j u^*_{0,j}e^{2 \pi i j n/N}$

### 5.1.6.1. Aliasing and nonlinearities¶

• assume we have a non-linear term $$\rho u$$ in our PDE, and
$\rho(z)=sin(k_1z),u(z)=sin(k_2z),$

with $$k_1,k_2$$ form our set of available wave-vectors $$k_j$$

• now
$\rho u ~ -cos[(k_1+k_2)z]+cos[(k_2-k_1)z],$

and $$k_1+k_2$$ may lie outside our range of k’s, and the available Foureier amplitudes might get aliased!

• $$k_1+k_2$$ outside range if $$k_1+k_2>\pi$$ and the amplitude appears wrongly in the range of ks at $$k_1+k_2-2\pi(l=-1,j_1+j_2-N)$$ the DFT is aliased

### 5.1.6.2. De-aliasing¶

• several methods exist to prevent aliasing:
• we apply 2/3-rule:
• simple to apply
• low cost in computing time
• Basic idea:

set part of the amplitudes to zero always prior to (non-linear) multiplications:

full index range of k-vectors: $$[-N/2,N/2]$$ -> keep the sub-range $$[-K,K]$$ free of aliasing method: set Fourier amplitudes $$u^*_j=0\ in \ [-N/2,-K]\ and \ [N/2, K]$$

• why does this work? and how to choose K?
• let j and s be in $$[0,K]$$
• if $$j+s>N/2$$ (outside range), then the amplitude corresponding to j+s will be aliased to j+s-N
• we demand that $$j+s-N<-K$$ (in the not used part of the spectrum), the largest j, s in the range are $$j=s=K; j+s-N<=2K-N$$ i.e. $$K<N/3$$
• we set $$K=N/3=(2/3)N/2$$: ‘2/3-rule`
• for $$j.s$$ in $$j+s>N/2$$ the amplitude is aliased to $$j+s-N$$, which may lie in $$[-K,0]$$ , but we do not have to care, the amplitudes at j and s are set to zero-> the range $$[-K,K]$$ is free of aliasing.

### 5.1.6.3. non-linearities, de-aliased¶

• assume you need to evaluate $$DFT(\rho_iu_i)$$ , having given the Fourier transforms $$\rho_j^*$$ and $$u_j^*$$ :
$\begin{split}&\rho_j^*,u_j^* \\ \rightarrow & \rho_j^*,u_j^*=0\ for\ j>(2/3)N/2 \\ \rightarrow &(DFT^{-1}) \rho_n , u_n \\ \rightarrow &w_n=\rho_n u_n \\ \rightarrow &(\rho_ju_j)^*=w^*\end{split}$

### 5.1.6.4. Stability and convergence¶

• ... theory on stability on convergence ...
• reproduce analytically known cases
• reproduce results of others, or results derived in different ways
• test the individual sub-tasks the code performs
• monitor conserved quantities (if there are any)
• apply fantasy and pysical intuition to the concrete problem you study, try to be as critical as you can against your results

## 5.1.7. Concluding remarks¶

### 5.1.7.1. Positive properties of the pseudo-spectral(PS) method¶

• for analytic functions (solutions), the errors decay exponentially with N, i.e. very fast
• non-smoothness of even discontinuities in coeficients ro the solutions seem note to cause problems
• often, less grid points are needed with the PS method than with finite difference methods to achieve the same accuracy(computing time and memory!)

### 5.1.7.2. Negative properties of the pseudo-spectral method¶

• certain boundary condition s may cause difficulties
• irregular domains(simulation boxes) can be difficult or impossible to implement
• strong shocks can cause problems
• local grid refinement (for cases where it is needed) seems not possible, so-far

## 5.1.8. References¶

 [1] Here is a nice tutorial
 [3] Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). “Section 20.7. Spectral Methods”. Numerical Recipes: The Art of Scientific Computing (3rd ed.) . New York: Cambridge University Press. ISBN 978-0-521-88068-8.