Spectral method




Spectral methods are a class of techniques used in applied mathematics and scientific computing to numerically solve certain differential equations, potentially involving the use of the fast Fourier transform. The idea is to write the solution of the differential equation as a sum of certain "basis functions" (for example, as a Fourier series which is a sum of sinusoids) and then to choose the coefficients in the sum in order to satisfy the differential equation as well as possible.


Spectral methods and finite element methods are closely related and built on the same ideas; the main difference between them is that spectral methods use basis functions that are nonzero over the whole domain, while finite element methods use basis functions that are nonzero only on small subdomains. In other words, spectral methods take on a global approach while finite element methods use a local approach. Partially for this reason, spectral methods have excellent error properties, with the so-called "exponential convergence" being the fastest possible, when the solution is smooth. However, there are no known three-dimensional single domain spectral shock capturing results (shock waves are not smooth).[1] In the finite element community, a method where the degree of the elements is very high or increases as the grid parameter h decreases to zero is sometimes called a spectral element method.


Spectral methods can be used to solve ordinary differential equations (ODEs), partial differential equations (PDEs) and eigenvalue problems involving differential equations. When applying spectral methods to time-dependent PDEs, the solution is typically written as a sum of basis functions with time-dependent coefficients; substituting this in the PDE yields a system of ODEs in the coefficients which can be solved using any numerical method for ODEs. Eigenvalue problems for ODEs are similarly converted to matrix eigenvalue problems[citation needed].


Spectral methods were developed in a long series of papers by Steven Orszag starting in 1969 including, but not limited to, Fourier series methods for periodic geometry problems, polynomial spectral methods for finite and unbounded geometry problems, pseudospectral methods for highly nonlinear problems, and spectral iteration methods for fast solution of steady state problems. The implementation of the spectral method is normally accomplished either with collocation or a Galerkin or a Tau approach.


Spectral methods are computationally less expensive than finite element methods, but become less accurate for problems with complex geometries and discontinuous coefficients. This increase in error is a consequence of the Gibbs phenomenon.




Contents





  • 1 Examples of spectral methods

    • 1.1 A concrete, linear example

      • 1.1.1 Algorithm



    • 1.2 Nonlinear example



  • 2 A relationship with the spectral element method


  • 3 See also


  • 4 References




Examples of spectral methods



A concrete, linear example


Here we presume an understanding of basic multivariate calculus and Fourier series. If g(x,y) is a known, complex-valued function of two real variables, and g is periodic in x and y (that is, g(x,y)=g(x+2π,y)=g(x,y+2π)) then we are interested in finding a function f(x,y) so that


(∂2∂x2+∂2∂y2)f(x,y)=g(x,y)for all x,ydisplaystyle left(frac partial ^2partial x^2+frac partial ^2partial y^2right)f(x,y)=g(x,y)quad textfor all x,yleft(fracpartial^2partial x^2+fracpartial^2partial y^2right)f(x,y)=g(x,y)quad textfor all x,y

where the expression on the left denotes the second partial derivatives of f in x and y, respectively. This is the Poisson equation, and can be physically interpreted as some sort of heat conduction problem, or a problem in potential theory, among other possibilities.


If we write f and g in Fourier series:


f=:∑aj,kei(jx+ky)displaystyle f=:sum a_j,ke^i(jx+ky)f=:sum a_j,ke^i(jx+ky)

g=:∑bj,kei(jx+ky)displaystyle g=:sum b_j,ke^i(jx+ky)g=:sum b_j,ke^i(jx+ky)

and substitute into the differential equation, we obtain this equation:


∑−aj,k(j2+k2)ei(jx+ky)=∑bj,kei(jx+ky)displaystyle sum -a_j,k(j^2+k^2)e^i(jx+ky)=sum b_j,ke^i(jx+ky)sum -a_j,k(j^2+k^2)e^i(jx+ky)=sum b_j,ke^i(jx+ky)

We have exchanged partial differentiation with an infinite sum, which is legitimate if we assume for instance that f has a continuous second derivative. By the uniqueness theorem for Fourier expansions, we must then equate the Fourier coefficients term by term, giving


(*) aj,k=−bj,kj2+k2displaystyle a_j,k=-frac b_j,kj^2+k^2a_j,k=-fracb_j,kj^2+k^2

which is an explicit formula for the Fourier coefficients aj,k.


With periodic boundary conditions, the Poisson equation possesses a solution only if b0,0 = 0. Therefore,
we can freely choose a0,0 which will be equal to the mean of the resolution. This corresponds to choosing the
integration constant.


To turn this into an algorithm, only finitely many frequencies are solved for. This introduces an error which can be shown to be proportional to hndisplaystyle h^nh^n, where h:=1/ndisplaystyle h:=1/nh:=1/n and ndisplaystyle nn is the highest frequency treated.



Algorithm


  1. Compute the Fourier transform (bj,k) of g.

  2. Compute the Fourier transform (aj,k) of f via the formula (*).

  3. Compute f by taking an inverse Fourier transform of (aj,k).

Since we're only interested in a finite window of frequencies (of size n, say) this can be done using a Fast Fourier Transform algorithm. Therefore, globally the algorithm runs in time O(n log n).



Nonlinear example


We wish to solve the forced, transient, nonlinear Burgers' equation using a spectral approach.


Given u(x,0)displaystyle u(x,0)u(x,0) on the periodic domain
x∈[0,2π)displaystyle xin left[0,2pi right)xinleft[0,2piright), find u∈Udisplaystyle uin mathcal Uu in mathcalU such that


∂tu+u∂xu=ρ∂xxu+f∀x∈[0,2π),∀t>0displaystyle partial _tu+upartial _xu=rho partial _xxu+fquad forall xin left[0,2pi right),forall t>0partial_t u + u partial_x u = rho partial_xx u + f quad forall xinleft[0,2piright), forall t>0

where ρ is the viscosity coefficient. In weak conservative form this becomes


⟨∂tu,v⟩=⟨∂x(−12u2+ρ∂xu),v⟩+⟨f,v⟩∀v∈V,∀t>0displaystyle leftlangle partial _tu,vrightrangle =leftlangle partial _xleft(-frac 12u^2+rho partial _xuright),vrightrangle +leftlangle f,vrightrangle quad forall vin mathcal V,forall t>0displaystyle leftlangle partial _tu,vrightrangle =leftlangle partial _xleft(-frac 12u^2+rho partial _xuright),vrightrangle +leftlangle f,vrightrangle quad forall vin mathcal V,forall t>0

where ⟨f,g⟩:=∫02πf(x)g(x)¯dxdisplaystyle langle f,grangle :=int _0^2pi f(x)overline g(x),dxlangle f, g rangle := int_0^2pi f(x) overlineg(x),dx"/> following inner product notation. Integrating by parts and using periodicity grants


⟨∂tu,v⟩=⟨12u2−ρ∂xu,∂xv⟩+⟨f,v⟩∀v∈V,∀t>0.displaystyle langle partial _tu,vrangle =leftlangle frac 12u^2-rho partial _xu,partial _xvrightrangle +leftlangle f,vrightrangle quad forall vin mathcal V,forall t>0.displaystyle langle partial _tu,vrangle =leftlangle frac 12u^2-rho partial _xu,partial _xvrightrangle +leftlangle f,vrightrangle quad forall vin mathcal V,forall t>0.

To apply the Fourier-Galerkin method, choose both


UN:=u:u(x,t)=∑k=−N/2N/2−1u^k(t)eikxdisplaystyle mathcal U^N:=leftu:u(x,t)=sum _k=-N/2^N/2-1hat u_k(t)e^ikxrightmathcalU^N := left u : u(x,t)=sum_k=-N/2^N/2-1 hatu_k(t) e^i k xright

and


VN:= spaneikx:k∈−N/2,…,N/2−1displaystyle mathcal V^N:=text spanlefte^ikx:kin -N/2,dots ,N/2-1rightmathcalV^N :=text spanleft e^i k x : kin -N/2,dots,N/2-1right

where u^k(t):=12π⟨u(x,t),eikx⟩displaystyle hat u_k(t):=frac 12pi langle u(x,t),e^ikxrangle hatu_k(t):=frac12pilangle u(x,t), e^i k x rangle. This reduces the problem to finding u∈UNdisplaystyle uin mathcal U^NuinmathcalU^N such that


⟨∂tu,eikx⟩=⟨12u2−ρ∂xu,∂xeikx⟩+⟨f,eikx⟩∀k∈−N/2,…,N/2−1,∀t>0.displaystyle langle partial _tu,e^ikxrangle =leftlangle frac 12u^2-rho partial _xu,partial _xe^ikxrightrangle +leftlangle f,e^ikxrightrangle quad forall kin left-N/2,dots ,N/2-1right,forall t>0.displaystyle langle partial _tu,e^ikxrangle =leftlangle frac 12u^2-rho partial _xu,partial _xe^ikxrightrangle +leftlangle f,e^ikxrightrangle quad forall kin left-N/2,dots ,N/2-1right,forall t>0.

Using the orthogonality relation ⟨eilx,eikx⟩=2πδlkdisplaystyle langle e^ilx,e^ikxrangle =2pi delta _lklangle e^i l x, e^i k x rangle = 2 pi delta_lk where δlkdisplaystyle delta _lkdelta_lk is the Kronecker delta, we simplify the above three terms for each kdisplaystyle kk to see


⟨∂tu,eikx⟩=⟨∂t∑lu^leilx,eikx⟩=⟨∑l∂tu^leilx,eikx⟩=2π∂tu^k,⟨f,eikx⟩=⟨∑lf^leilx,eikx⟩=2πf^k, and⟨12u2−ρ∂xu,∂xeikx⟩=⟨12(∑pu^peipx)(∑qu^qeiqx)−ρ∂x∑lu^leilx,∂xeikx⟩=⟨12∑p∑qu^pu^qei(p+q)x,ikeikx⟩−⟨ρi∑llu^leilx,ikeikx⟩=−ik2⟨∑p∑qu^pu^qei(p+q)x,eikx⟩−ρk⟨∑llu^leilx,eikx⟩=−iπk∑p+q=ku^pu^q−2πρk2u^k.displaystyle beginalignedleftlangle partial _tu,e^ikxrightrangle &=leftlangle partial _tsum _lhat u_le^ilx,e^ikxrightrangle =leftlangle sum _lpartial _that u_le^ilx,e^ikxrightrangle =2pi partial _that u_k,\leftlangle f,e^ikxrightrangle &=leftlangle sum _lhat f_le^ilx,e^ikxrightrangle =2pi hat f_k,text and\leftlangle frac 12u^2-rho partial _xu,partial _xe^ikxrightrangle &=leftlangle frac 12left(sum _phat u_pe^ipxright)left(sum _qhat u_qe^iqxright)-rho partial _xsum _lhat u_le^ilx,partial _xe^ikxrightrangle \&=leftlangle frac 12sum _psum _qhat u_phat u_qe^ileft(p+qright)x,ike^ikxrightrangle -leftlangle rho isum _llhat u_le^ilx,ike^ikxrightrangle \&=-frac ik2leftlangle sum _psum _qhat u_phat u_qe^ileft(p+qright)x,e^ikxrightrangle -rho kleftlangle sum _llhat u_le^ilx,e^ikxrightrangle \&=-ipi ksum _p+q=khat u_phat u_q-2pi rho k^2hat u_k.endaligneddisplaystyle beginalignedleftlangle partial _tu,e^ikxrightrangle &=leftlangle partial _tsum _lhat u_le^ilx,e^ikxrightrangle =leftlangle sum _lpartial _that u_le^ilx,e^ikxrightrangle =2pi partial _that u_k,\leftlangle f,e^ikxrightrangle &=leftlangle sum _lhat f_le^ilx,e^ikxrightrangle =2pi hat f_k,text and\leftlangle frac 12u^2-rho partial _xu,partial _xe^ikxrightrangle &=leftlangle frac 12left(sum _phat u_pe^ipxright)left(sum _qhat u_qe^iqxright)-rho partial _xsum _lhat u_le^ilx,partial _xe^ikxrightrangle \&=leftlangle frac 12sum _psum _qhat u_phat u_qe^ileft(p+qright)x,ike^ikxrightrangle -leftlangle rho isum _llhat u_le^ilx,ike^ikxrightrangle \&=-frac ik2leftlangle sum _psum _qhat u_phat u_qe^ileft(p+qright)x,e^ikxrightrangle -rho kleftlangle sum _llhat u_le^ilx,e^ikxrightrangle \&=-ipi ksum _p+q=khat u_phat u_q-2pi rho k^2hat u_k.endaligned

Assemble the three terms for each kdisplaystyle kk to obtain


2π∂tu^k=−iπk∑p+q=ku^pu^q−2πρk2u^k+2πf^kk∈−N/2,…,N/2−1,∀t>0.displaystyle 2pi partial _that u_k=-ipi ksum _p+q=khat u_phat u_q-2pi rho k^2hat u_k+2pi hat f_kquad kin left-N/2,dots ,N/2-1right,forall t>0.<br/>2 pi partial_t hatu_k<br/>=<br/>- i pi k sum_p+q=- 2pirhok^2hatu_k
+ 2 pi hatf_k
quad kinleft -N/2,dots,N/2-1 right, forall t>0.
"/>

Dividing through by 2πdisplaystyle 2pi 2pi , we finally arrive at


∂tu^k=−ik2∑p+q=ku^pu^q−ρk2u^k+f^kk∈−N/2,…,N/2−1,∀t>0.displaystyle partial _that u_k=-frac ik2sum _p+q=khat u_phat u_q-rho k^2hat u_k+hat f_kquad kin left-N/2,dots ,N/2-1right,forall t>0.<br/>partial_t hatu_k<br/>=<br/>- fraci k2 sum_p+q=- rhok^2hatu_k
+ hatf_k
quad kinleft -N/2,dots,N/2-1 right, forall t>0.
"/>

With Fourier transformed initial conditions u^k(0)displaystyle hat u_k(0)hatu_k(0) and forcing f^k(t)displaystyle hat f_k(t)hatf_k(t), this coupled system of ordinary differential equations may be integrated in time (using, e.g., a Runge Kutta technique) to find a solution. The nonlinear term is a convolution, and there are several transform-based techniques for evaluating it efficiently. See the references by Boyd and Canuto et al. for more details.



A relationship with the spectral element method


One can show that if gdisplaystyle gg is infinitely differentiable, then the numerical algorithm using Fast Fourier Transforms will converge faster than any polynomial in the grid size h. That is, for any n>0, there is a Cn<∞displaystyle C_n<infty displaystyle C_n<infty such that the error is less than Cnhndisplaystyle C_nh^ndisplaystyle C_nh^n for all sufficiently small values of hdisplaystyle hh. We say that the spectral method is of order ndisplaystyle nn, for every n>0.


Because a spectral element method is a finite element method of very high order, there is a similarity in the convergence properties. However, whereas the spectral method is based on the eigendecomposition of the particular boundary value problem, the finite element method does not use that information and works for arbitrary elliptic boundary value problems.



See also


  • Finite element method

  • Gaussian grid

  • Pseudo-spectral method

  • Spectral element method

  • Galerkin method

  • Collocation method



References




  1. ^ pp 235, Spectral Methods: evolution to complex geometries and applications to fluid dynamics, By Canuto, Hussaini, Quarteroni and Zang, Springer, 2007.



  • Bengt Fornberg (1996) A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge, UK


  • Chebyshev and Fourier Spectral Methods by John P. Boyd.

  • Canuto C., Hussaini M. Y., Quarteroni A., and Zang T.A. (2006) Spectral Methods. Fundamentals in Single Domains. Springer-Verlag, Berlin Heidelberg

  • Javier de Frutos, Julia Novo: A Spectral Element Method for the Navier--Stokes Equations with Improved Accuracy


  • Polynomial Approximation of Differential Equations, by Daniele Funaro, Lecture Notes in Physics, Volume 8, Springer-Verlag, Heidelberg 1992

  • D. Gottlieb and S. Orzag (1977) "Numerical Analysis of Spectral Methods : Theory and Applications", SIAM, Philadelphia, PA

  • J. Hesthaven, S. Gottlieb and D. Gottlieb (2007) "Spectral methods for time-dependent problems", Cambridge UP, Cambridge, UK

  • Steven A. Orszag (1969) Numerical Methods for the Simulation of Turbulence, Phys. Fluids Supp. II, 12, 250-257


  • 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..mw-parser-output cite.citationfont-style:inherit.mw-parser-output .citation qquotes:"""""""'""'".mw-parser-output .citation .cs1-lock-free abackground:url("//upload.wikimedia.org/wikipedia/commons/thumb/6/65/Lock-green.svg/9px-Lock-green.svg.png")no-repeat;background-position:right .1em center.mw-parser-output .citation .cs1-lock-limited a,.mw-parser-output .citation .cs1-lock-registration abackground:url("//upload.wikimedia.org/wikipedia/commons/thumb/d/d6/Lock-gray-alt-2.svg/9px-Lock-gray-alt-2.svg.png")no-repeat;background-position:right .1em center.mw-parser-output .citation .cs1-lock-subscription abackground:url("//upload.wikimedia.org/wikipedia/commons/thumb/a/aa/Lock-red-alt-2.svg/9px-Lock-red-alt-2.svg.png")no-repeat;background-position:right .1em center.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registrationcolor:#555.mw-parser-output .cs1-subscription span,.mw-parser-output .cs1-registration spanborder-bottom:1px dotted;cursor:help.mw-parser-output .cs1-ws-icon abackground:url("//upload.wikimedia.org/wikipedia/commons/thumb/4/4c/Wikisource-logo.svg/12px-Wikisource-logo.svg.png")no-repeat;background-position:right .1em center.mw-parser-output code.cs1-codecolor:inherit;background:inherit;border:inherit;padding:inherit.mw-parser-output .cs1-hidden-errordisplay:none;font-size:100%.mw-parser-output .cs1-visible-errorfont-size:100%.mw-parser-output .cs1-maintdisplay:none;color:#33aa33;margin-left:0.3em.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registration,.mw-parser-output .cs1-formatfont-size:95%.mw-parser-output .cs1-kern-left,.mw-parser-output .cs1-kern-wl-leftpadding-left:0.2em.mw-parser-output .cs1-kern-right,.mw-parser-output .cs1-kern-wl-rightpadding-right:0.2em

  • Lloyd N. Trefethen (2000) Spectral Methods in MATLAB. SIAM, Philadelphia, PA


Popular posts from this blog

How to check contact read email or not when send email to Individual?

Displaying single band from multi-band raster using QGIS

How many registers does an x86_64 CPU actually have?