3.1. Pseudo Spectral Method 算法说明

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

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

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)

3.1.3. Pseudo-spectral methods and Fourier transforms 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:
  • 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\]


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


(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

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

3.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\)

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

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^*\) :
\[\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}\] 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

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

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