4.1. Pseudo Spectral Method 算法说明

4.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
  • ...

4.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:

L u(x)=s(x), x\in V \\
b.c.:f(u(y))=0, y\in \partial V

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

  • wanted: numerical solution u^N(x) such that the residual R:


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


  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)

4.1.3. Pseudo-spectral methods and Fourier transforms Comparison to analytical Fourier method

  • in

Lu(x)=s(x),\ \ x\in V\\
b.c.:f(u(y))=0,\ y\in \partial V

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

  • Fourier transform:


  • in principle, we want to do this numerically, but we have to make sure about a few points... 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


\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


k_j= 2\pi j/(N \Delta), j=-N/2,\dots,N/2\\
\Delta:\textrm{spatial resotuion}

  • 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


(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

  • 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

4.1.4. Time-dependent problems

  • example: diffusion equation in 1D:


  • we consider the euqation only at the collocation points z_n=n\Delta,n=0,1,\dots,N-1 writing symbolically


apply a spatial DFT:


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

we have a set of N ODEs the temporal integration is done in Fourier space 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:

du/dt&=F(u,t) \ \ \ (u has N components)\\
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)\\

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

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

4.1.6. Aliasing–de-aliansing

  • the Foureier modes used are:

\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}

  • 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


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} Aliasing and nonlinearities

  • assume we have a non-linear term \rho u in our PDE, and


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 De-aliasing

  • several methods exist to prevent aliasing:
zero-padding(3/2-rule), truncating(2/3-rule), phase shift
  • 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. non-linearities, de-aliased

  • assume you need to evaluate DFT(\rho_iu_i) , having given the Fourier transforms \rho_j^* and u_j^* :

&\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^* 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

4.1.7. Concluding remarks 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!) 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

4.1.8. References

[1]Here is a nice tutorial
[2]Wikipedia page for pseudo spectral method
[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.