Spectral method
|
|
In applied mathematics, spectral methods are algorithms to solve certain kinds of partial differential equations numerically using some sort of Fast Fourier Transform.
Partial differential equations (PDEs) describe a wide array of physical processes such has heat conduction and sound propagation. In many such equations, there are underlying "basic waves" that can be used to give efficient algorithms for computing solutions to these PDEs.
In all cases, one is provided with input data (such as a description of the medium and initial condition for a sound propagation experiment) and a specific PDE governing the physical process. One then writes the input data as a superposition of the basic waves of the PDE. In this form, the solution of the PDE is easily calculated.
See also: Finite element method Discrete element method
A concrete example
Here we presume a basic understanding of basic multivariate calculus and Fourier series. If g(x,y) is a 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
- <math>\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)f(x,y)=g(x,y)\quad for\ all\ x,y<math>
where the things on the left denote the second partial derivatives of f in x and y, respectively. This is Laplace's equation, and can be physically interpreted as some sort of heat conduction problem.
If we write f and g in Fourier series:
- <math>f=\sum a_{j,k}e^{ijx+iky}<math>
- <math>g=\sum b_{j,k}e^{ijx+iky}<math>
and substitute into the differential equation, we obtain this equation:
- <math>-a_{j,k}(j^2+k^2)e^{ijx+iky}=\sum b_{j,k}e^{ijx+iky}<math>
Here, we crashed the differentials through the infinite sums, but a justification is beyond the scope of this article. By the uniqueness theorem for Fourier expansions, we must then equate the Fourier coefficients term by term, giving
- (*) <math>a_{j,k}=-\frac{b_{j,k}}{j^2+k^2}<math>
which is an explicit formula for the Fourier coefficients aj,k.
To turn this into an algorithm, only finitely many frequencies are solved for.
Algorithm
- Compute the Fourier transform (bj,k) of g.
- Compute the Fourier transform (aj,k) of f via the formula (*) and the Fourier transform of g.
- 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(nlogn).