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)\]

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

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.

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

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