While the general strategy of Partial Differential Equation Toolbox™ software is to use the MATLAB® built-in solvers for sparse systems, there are situations where faster solution algorithms are available. One such example is found when solving Poisson's equation
–Δu = f in Ω
with Dirichlet boundary conditions, where Ω is a rectangle.
For the fast solution algorithms to work, the mesh on the rectangle must be a regular mesh. In this context it means that the first side of the rectangle is divided into N1 segments of length h1, the second into N2 segments of length h2, and (N1 + 1) by (N2 + 1) points are introduced on the regular grid thus defined. The triangles are all congruent with sides h1, h2 and a right angle in between.
The Dirichlet boundary conditions are eliminated in the usual way, and the resulting problem for the interior nodes is Kv = F. If the interior nodes are numbered from left to right, and then from bottom to top, the K matrix is block tridiagonal. The N2 – 1 diagonal blocks, here called T, are themselves tridiagonal (N1 – 1) by (N1 – 1) matrices, with 2(h1/h2 + h2/h1) on the diagonal and –h2/h1 on the subdiagonals. The subdiagonal blocks, here called I, are –h1/h2 times the unit N1 – 1 matrix.
The key to the solution of the problem Kv = F is that the problem Tw = f is possible to solve using the discrete sine transform. Let S be the (N1 – 1) by (N1 – 1) matrix with Sij= sin(πij/N1). Then S–1TS = Λ, where Λ is a diagonal matrix with diagonal entries 2(h1/h2 + h2/h1) – 2h2/h1 cos(πi/N1). w = SΛ–1S–1 f, but multiplying with S is nothing more than taking the discrete sine transform, and multiplying with S–1 is the same as taking the inverse discrete sine transform. The discrete sine transform can be efficiently calculated using the fast Fourier transform on a sequence of length 2N1.
Solving Tw = f using the discrete sine transform would not be an advantage in itself, since the system is tridiagonal and should be solved as such. However, for the full system Ky = F, a transformation of the blocks in K turns it into N1 – 1 decoupled tridiagonal systems of size N2 – 1. Thus, a solution algorithm would look like
When using a fast solver such as this one, time and memory are also saved since the matrix K in fact never has to be assembled. A drawback is that since the mesh has to be regular, it is impossible to do adaptive mesh refinement.
The fast elliptic solver for Poisson's equation is implemented in poisolv. The discrete sine transform and the inverse discrete sine transform are computed by dst and idst, respectively.